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ABSTRACT 



Passive microwave satellite data provides all-weather, continuous 
observations of the polar ice cover. Multichannel radiometer data sets extending 
back to 1978 are readily available and distributed on CD-ROM media. Current 
algorithms can extract information on ice cover, concentration and type. Two 
applications of these data sets for polar studies are explored: 1) the possible 

detection of altered ice morphology due to anomalous tidal currents in the 
Yermak Plateau region; and 2) the tracking of multiyear ice features to infer ice 
motion. 

No evidence of altered ice morphology in the Yermak Plateau region was 
detected. Either the reported anomalous currents do not affect ice morphology or 
the limitations inherent to passive microwave imagery (poor spatial resolution) 
prevent detection. On the other hand, distinct multiyear ice features are observed 
to move in the direction of mean ice flow and are well correlated with the motion 
of nearby data buoys. Short term ice motion vectors are unreliable because of the 
low spatial resolution of the passive microwave data. However, long term 
averages are indicative of ice flow trajectories and can be used, for example, to 
increase our knowledge of the interannual variability of heat and salinity budgets 
in specific regions, including Fram Strait, the principal outflow region of the arctic 
basin. 
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I. INTRODUCTION 



All-weather, continuous observations of the arctic ice cover are important for 
commercial and naval operations and for scientific studies such as global climate 
change. However, such observations are generally difficult to obtain. In-situ 
measurements from ice camps and ships are expensive to acquire and spatially 
limited. Aircraft measurements can expand the areal coverage, but asset and 
resource limitations prevent continuous large-scale monitoring of the ice pack. 
Satellite platforms are the only practical means of providing continuous coverage. 
Of the passive sensing methods, only microwave emissions penetrate the clouds 
and darkness that shield much of the arctic ice from visual and infrared 
observations. Passive satellite microwave data are operationally available from 
Special Sensor Microwave/Imager (SSM/I) sensors on the Defense 
Meteorological Satellite Program (DMSP) series satellites (Hollinger et ah, 1987). 
Historical data are also available to researchers on CD-ROM distributed by the 
National Snow and Ice Data Center (NSIDC). 

The objectives of this study are to explore the use of passive microwave 
satellite measurements for arctic sea ice studies. Two applications are examined: 

• The detection of altered ice morphology due to anomalous 
diurnal tidal currents observed over the Yermak Plateau. If 
successful, this technique can be applied in other polar regions to 
expand our knowledge of the local tidal structure and its effect on 
sea ice cover. 

• Determination of ice flow by tracking the movement of floes 
exhibiting different morphologies from the surrounding ice field. 
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This can augment the Arctic Buoy Program and increase our 
understanding of ice dynamics. 

Chapter II presents background material on the SSM/I, data sources and sea 
ice microwave radiation characteristics. Determination of sea ice concentrations 
from passive microwave data are also discussed. Chapter III examines the 
anomalous diurnal currents observed over the Yermak Plateau. Previous 
observations and models prove the existence and explain the physics behind 
these currents. The effect of these currents on the ice cover is examined through 
passive microwave emissions. Chapter IV shows how ice features can be tracked 
to infer general ice drift patterns in the arctic. When combined with other data, 
this technique can help researchers understand ice and salinity budgets in the 
polar region. Chapter V summarizes the results and closes with a discussion of 
further research that could be attempted using these data. 
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II. BACKGROUND 



A. THE SSM/I 

The SSM/I is flown on the Defense Meteorological Satellite Program (DMSP) 
block 5D-2 F8 satellite, launched on June 19, 1987. It is in circular sun- 
synchronous near-polar orbit at an altitude of approximately 833 km, an 
inclination of 98.8° and a period of 102 minutes. This produces 14.1 full orbit 
revolutions per day with a 0613 local ascending node equatorial crossing. The 
SSM/I's conical scan covers a 1400 km swath which provides 24 hour near-global 
coverage. The orbit inclination precludes coverage of a 280 km radius region 
centered on each pole extending to about 87.5°N. Equatorial regions require 72 
hours for full coverage, however overlapping orbits provide full daily coverage 
poleward of 45°N. Hollinger et al. (1990) provide a good discussion of SSM/I 
operating characteristics. 

The SSM/I uses four frequencies to measure seven channels of linearly 
polarized radiation as shown in Table 1. Footprint resolution increases with 
increasing frequency. Unfortunately, the high resolution 85 GHz vertically 
polarized (85V) channel was damaged when the satellite overheated in the winter 
of 1987. The channel was useless by the end of January 1989 (Hollinger et al., 
1990). The National Snow and Ice Data Center (NSIDC) recently reported that 
the 85 GHz horizontally polarized (85H) channel is becoming increasingly 
unreliable and is basically useless (Taylor, 1990). Poe and Conway (1990) 
reported the current SSM/I geo-location error to be approximately 10-15 km. 
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NSIDC applies a correction algorithm that reduces this error to 8 km (NSIDC, 
1990). 



TABLE 1: SSM/I FREQUENCIES 



Channel 


Frequency 

(GHz) 


Wavelength 

(cm) 


Polarization 


Resolution 

(km) 


1 


19.35 


1.55 


Vertical 


56 


2 


19.35 


1.55 


Horizontal 


56 


3 


22.235 


1.35 


Vertical 


45 


4 


37.0 


0.81 


Vertical 


33 


5 


37.0 


0.81 


Horizontal 


33 


6 


85.5 


0.35 


Vertical 


14 


7 


85.5 


0.35 


Horizontal 


14 



SSM/I data are transmitted to Telemetry and Data Relay Stations on each orbit. 
These data are then relayed via satellite to Navy and Air Force processing centers 
which produce near real time global maps of sea ice location, age, type and other 
environmental parameters. Raw data are also continually broadcast to Air Force 
and Navy tactical units positioned within the orbital swath of the satellite that are 
equipped with SMQ-11 satellite receivers. 
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Data record files are created from raw radiometer measurements. Temperature 
data records (TDRs) are antenna temperatures which are calibrated each scan 
using a 300° K black-body radiator and the 3° K cosmic radiation as reference 
temperatures. An antenna pattern correction is applied to the TDRs to correct for 
energy received in the sidelobes, cross-polarization coupling and feedhorn 
spillover loss (NSIDC, 1990). These main beam brightness temperatures are 
placed in sensor data record (SDR) files. Environmental parameters are 
determined by applying appropriate algorithms to the SDRs and are stored in 
environmental data records (EDRs). The primary archiver of SSM/I data in orbital 
format is the NOAA/National Environmental Satellite, Data, and Information 
Service, Satellite Data Services Division (NOAA/NESDIS/SDSD). TDRs, SDRs 
and EDRs are also forwarded to the Naval Research Laboratory for archival. 

B. SEA ICE RADIATION CHARACTERISTICS 

Sea ice is a complex material. A good comprehensive overview of both its 
physical and microwave radiation characteristics can be found in Parkinson et al. 
(1987). The following is a general discussion of these properties. 

1. Physical Properties 

Sea ice covers approximately seven percent of the word's oceans. It is 
inhomogeneous, anisotropic and never in equilibrium. Sea ice can simultaneously 
contain all three phases of matter: solid ice, liquid brine pockets and air cavities. 
The fractions of the constituent phases change over time. New ice is thinner and 
saltier than first year ice, which transforms after a melt cycle (summer season) into 
fresher, thicker and less dense multiyear ice containing numerous air pockets. Ice 
can be covered with a layer of dry or wet snow of varying thickness; summer 
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melting forms ponds of liquid water over the ice sheet. Figure 1 shows the 
different physical characteristics of first and multiyear ice. 

The ice composition is initially determined by environmental factors at 
the time it is formed. These factors include air temperature, wind speed, ocean 
currents, near surface salinity and weather (Weeks, 1981). The ice composition 
changes continuously due to weather, ice advection and fracturing or ridging 
effects. Warmer summer temperatures are especially important as melt water 
percolates through the ice and brine pockets are flushed. 

2. Radiative Properties 

The radiative properties of sea ice are determined by the emissivity 
differences between ice and water. Emissivity is the ratio of emittance (J nr 2 s" 1 ) 
per unit wavelength to that of a blackbody at the same temperature. The 
Rayleigh-Jeans approximation (linear relationship between emittance and 
temperature) applies at microwave wavelengths. The emissivity (e) of a substance 
can be expressed as: 



£ = 



T e 

Tsfc 



( 1 ) 



where T e is the radiated microwave brightness temperature, proportional to the 
emittance, and T s fc is the physical temperature of the substance. The SSM/I 
radiometer senses Tb, the brightness temperature at the satellite which is related 
to Te after attenuation and contamination by atmospheric effects are considered. 

Brightness temperatures can be used to calculate ice concentrations and 
differentiate between open water, first year ice and multiyear ice based on 
emissivity variations of the substances. To a lesser extent, snow and water- 
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covered ice can also be observed. Figure 2 summarizes the emissivity differences 
between sea water, first year ice and multiyear ice as a function of polarized 
frequency. 

The ice composition controls the dielectric properties that influence the 
emissivity of the ice. Following Comiso (1983), the principal difference between 
sea ice and freshwater ice are brine inclusions in the former. These brine pockets 
are small, typically 3-5 mm by 0.25 mm (Vant et al., 1978). The brine content of 
sea ice depends on freezing conditions and subsequent temperature changes and 
varies from 16 psu for first year ice to less than 1 psu for multiyear ice. The liquid 
brine effects the dielectric property that determines emissivity and the skin depth, 
which is the depth from which most of the observed radiation emanates. In saline 
ice the skin depth is on the order of the radiation wavelength. For desalinated ice 
this depth can be many wavelengths (Gloersen et al., 1973). 

The fundamental difference between first year and multiyear ice is the 
presence of brine pockets in the former and their replacement in the freeboard 
layer by air in the latter. Because of the small skin depth, microwave emissions 
from the more saline first year ice emanate from the snow/ice interface; emissions 
from fresher multiyear ice emanate from a thicker ice layer and are subject to 
scattering by the air pockets. This scattering reduces the net radiation from the 
ice surface thus lowering the brightness temperature (Comiso, 1983). 

Snow cover also effects microwave radiation. Snow insulates the ice 
from the atmosphere, leading to warmer ice from oceanic heat conduction and 
therefore slightly higher brightness temperatures. This may be compensated by 
the scattering effect of snow ice crystals, especially at higher frequencies. 
Emissions from the crystals are negligible (Comiso et al., 1982). During the 
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summer melt season, however, each crystal is surrounded by a thin layer of water. 
This wet snow cover becomes the primary radiation source, significantly lowering 
brightness temperatures due to water’s low emissivity. This effect makes it nearly 
impossible to discriminate wet snow-covered multiyear ice from first year ice in 
the summer. 

Ice thickness also influences brightness temperatures. At microwave 
frequencies, thin ice (radiation wavelength or less) is contaminated by the coldly 
radiating water below, lowering the brightness temperature. This may be 
compensated by the higher salinity and warmer physical temperature of the thin 
ice versus the relatively fresher and colder thick ice. 

Melt ponds are formed by solar heating of the ice cover and are 
characterized by thin layers of water scattered among solid ice regions. The size 
and depth of the ponds depend on the incident solar flux, ice salinity and ice 
roughness (Goroch and Fett, 1990) If the thin water layer is thicker than the skin 
depth of the radiation (typically a few mm to one cm), the radiation is emitted 
primarily by the liquid and has a much colder brightness temperature than the ice. 
For thinner layers the underlying ice will contribute some radiation, raising the 
brightness temperature. The variation of brightness temperature depends on the 
layer depth and the emission frequency. For depths on the order of a few 
millimeters the 19 and 22 GHz channels will be unaffected (Goroch and Fett, 
1990). As the frequency increases, contamination by the water layer also 
increases. The total effect on brightness temperature is a function of areal extent 
within the radiometer footprint and melt pond depth. 

Other factors affecting microwave emissivities include the physical 
temperature of ice and undulations caused by ridging and refreezing of melt 
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ponds. As with liquid melt ponds, these conditions vary spatially and temporally 
and could cause observable changes in ice emissivity. 

3. Atmospheric Effects 

Important atmospheric constituents which affect microwave radiation 
are water vapor, oxygen and liquid water content. Seasonal climatological 
constants developed from physical models or observations are generally used in 
algorithms to quantify their cumulative effect on brightness temperatures. 

Interaction with non-precipitating clouds in the arctic region is 
negligible due to the small size of the droplets compared with the microwave 
wavelengths (Comiso, 1983). Liquid precipitation, however, can contaminate 
brightness temperatures. Operational algorithms usually flag rain-contaminated 
radiances for rejection. 

Although the arctic regions have low humidities, absorption by water 
vapor can effect remotely sensed emissions. For the SSM/I only the 22 and 85.5 
GHz channels are significantly effected by water vapor absorption (Hollinger et 
al, 1987). Neither channel is currently used in ice studies, although the higher 
resolution offered by operational 85.5 GHz channels in future DMSP satellites 
will be tempting. 

The atmospheric oxygen distribution is well known and constant. 
Absorption effects are included in seasonal climatological constants used in ice 
concentration algorithms. 

C. SEA ICE CONCENTRATION ALGORITHM 
1 . Physical Model 

A simple radiative model expressing the radiation sensed at satellite 
height (Tg) is presented in Svendsen et al. (1983). Surface and atmospheric 
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effects are included, plus a minor free space contribution. The model is composed 
of four terms (see Figure 3): radiation from the surface as seen through the 

atmosphere (A); upwelled radiation from the atmosphere (B); downwelled 
atmospheric radiation reflected at the surface (C); and free space radiation 
reflected by the surface (D). 



Tr — e sfcTsfc(l' x a) + ^T a x a + ( 1 ' e sfc)^T a x a ( 1 " T a) + ^"^sfc^sp (^) 
(A) (B) (C) (D) 

x a is the optical depth, T S p is the free space temperature and 8T a is the weighted 

average atmospheric temperature in the lower troposphere. Since the optical 
depth is small, the extinction exponential (e _x a) is approximated by l-x a . 
Atmospheric absorption of free space radiation is neglected. 

From equation (1) the emitted brightness temperature (T e ) is dependent 
only on the physical temperature of the emitter (T s f c ) and its emissivity (e s fc)- By 
substituting Te/T s f c for e s f c in equation (2), T e can be computed if x a and 8T a are 
known, and if T s f c is inferred from 5T a . 

Three components contribute to T e : sea water, first year ice and 
multiyear ice. The sum of the fractional components of each (C w , Cf and C m) 
respectively) add up to 1 : 

1 = C m + Cf + C w (3) 

The emitted brightness temperature is the sum of the individual components, 
weighted by the concentrations: 



10 



e sfc^sfc ^'m £ m^ice + ^-'f € f^ice + 



(4) 



where water temperature is assumed to be 272° K and £f, £ m and £ w are the 
emissivities of first year ice, multiyear ice and seawater, respectively. Tj C e in 

equation (4) is determined by assuming a linear relationship between the 
atmospheric surface temperature (T a ) and the water temperature: 

T ice = aT a + (l-a)(272) (5) 

where a is an observationally determined constant 0.4 from Svendsen et al., 
1983). 

Fractional concentrations can be computed from equations (1) - (5) 
using two radiometer channels. Sources of error include approximations of 
optical depth, ice temperature, atmospheric temperature weighting factor, 
atmospheric surface temperatures and assuming constant emissivities for first and 
multiyear ice. Svendsen et al. (1983) compared aircraft radiometric measurements 
with Nimbus 7 radiometer measurements using this technique with in-situ 
temperature observations and determined accuracies to be ± 3% of total ice 
concentration and ± 10% for multiyear ice concentration with a spatial resolution 
of 90 km. 

2. Multichannel Algorithm 
a. Multispectral Analysis 

Comiso (1983) applied a multi-spectral technique that exploits the 
variation of emissivity with frequency for sea water and various ice types to 
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distinguish between various first and multiyear ice and open water. Figure 4 
shows three distinct brightness temperature clusters that, when geo-located, 
correspond to first year ice, a first/multiyear ice mixture and open ocean 
conditions. The line that connects the first year cluster with the open ocean is 
interpreted as varying ice concentrations. Because the skin depth is frequency 
dependent, the grouping beneath the first year ice is theorized to be thin ice, 
where cold contamination from sea water is observed in the lower (deeper 
penetrating) frequency. The larger spread of emissivities and brightness 
temperatures in the 37V channel reflects the greater sensitivity of 37 GHz 
radiation to near-surface effects due to its smaller skin depth. 
b . Polarization! Gradient Ratios 

Cavelieri et al. (1984) used polarization ratios (PR) for a single 
frequency and spectral gradient ratios (GR) between two frequencies to estimate 
ice concentrations and fractions. PR and GR are defined as follows: 



T B v(fl)-T B h(fl) 

-Tb#!)^^) 



( 6 ) 



T B y( f l) ~ T B y( f 2) 
~ T Bv (fl) + T B y( f 2) 



( 7 ) 



where T Bv and Tg^ are brightness temperatures for vertical and horizontal 
polarizations respectively and fj and f 2 are separate frequencies. Table 2 shows 

PRs and GRs for sea ice and water emissivities at the two SSM/I frequencies (37.0 
and 19.35 GHz) used to calculate ice concentrations and fractions. Note that PR 
is much greater over the open ocean than over sea ice. Thus PR is a good 
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water/ice discriminator and can be used to estimate ice concentrations. GR 
depends on ice type and is affected by season. Ice types are determined by the 
sign of GR: positive over the open ocean, negative over multiyear ice and near 
zero over first year ice (Parkinson and Cavalieri, 1989). 



TABLE 2: EMISSIVITIES AT 19.35 AND 37 GHZ 



Ice Type 


First Year 
Thin Ice 


First Year 
Ice 


Multiyear 

Ice 


Calm 

Sea Water 


Polarization 


V 


H 


V 


H 


V 


H 


V 


H 


Frequency 

(GHz) 


Hi: 

■; • • • . \ - v : . • < £ . • ; :y; .• 

: 


19.35 


0.96 


0.78 


0.97 


0.84 


0.86 


0.73 


0.60 


0.33 


37.0 


0.96 


0.87 


0.97 


0.95 


0.69 


0.64 


0.70 


0.40 


Polarization 

Ratio 

(37V - 37H) 


0.09 


0.02 


0.05 


0.30 


Gradient 

Ratio 

(37V- 19V) 


0.00 


0.00 


i 

O 

* 


0.10 



Cavalieri et al. (1984) describe how these ratios are used to 
calculate total and multiyear ice concentrations and the physical temperature of 
the sea surface. They conclude that the precision of the calculated ice 
concentrations ranges from 5 - 9 % and the precision of the calculated multiyear 
ice fractions vary from 13 - 25 %. A caveat on the ice fraction calculation is that 
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first and multiyear ice can not be discriminated during the summer season 
because the general surface melt conditions mask the distinguishing properties of 
the two ice types. A non-iterative version of this method has been adopted by 
the NASA Sea Ice Algorithm Working Group and is used by the NSIDC (NSIDC, 
1990). 

D. DATA SOURCES 
1. NSIDC CD-ROMs 

NSIDC distributes SSM/I brightness temperatures on CD-ROM disks for 
both arctic and antarctic regions. A continuous record is currently available from 
9 July 1987 - 31 December 1989. Earlier SMMR data are available in the same 
format from 1978-1987. Subsequent SSM/I data are in production for CD-ROM 
distribution. The data include all SSM/I channels with 85 GHz data gridded at 
12.5 km resolution and 19, 22 and 37 GHz data gridded at 25 km resolution. 
Added features are land masks and coast outlines that provide geographic 
reference points on the imagery. The data are obtained from 
NOAA/NESDIS/SDSD. Fleet Numerical Oceanography Center produces the 
original TDR files from the raw satellite data as part of their operational data 
processing. The TDRs are transferred to NOAA/NESDIS/SDSD and then to 
NSIDC. At NSIDC the TDRs are examined for bad data values and a geo- 
location correction is applied. The antenna temperatures are then converted to 
brightness temperatures. 

Swath data from each channel are mapped onto the appropriate grid by 
converting the geodetic latitudes and longitudes at the center of each 
observation footprint into map grid coordinates. Brightness temperatures are 
summed in each grid cell over a 24 hour period and then divided by the number 
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of observations. The resulting data fields are daily averaged mosaics of each 
polar region. 

The data are mapped onto polar stereographic projections with the grid 
tangent to the earth at 70° latitude rather than at the poles. This increases the 
distortion at the poles by 3 % but is offset by nearly true projections in the 
marginal ice zones. Figure 5 shows the area of coverage for the Northern 
Hemisphere. 

2. Ice Buoy Tracks 

The Arctic Buoy Program has been operating since 1979. Objectives of 
the program are to provide surface atmospheric pressure measurements over the 
arctic basin and to define the large scale field of sea ice motion (Thorndike and 
Colony, 1980). Buoys are parachuted from C-130 aircraft at designated 
locations. Each buoy transmits pressure and air temperature information to a 
satellite at one minute intervals. This information is only received when the 
satellite is within the radio field of view. The data are retransmitted to a ground 
station, decoded and filtered, and forwarded to the Polar Science Center at the 
University of Washington. Thorndike and Colony (1980) examined the 
locational record of one buoy and reported the positional accuracy, as measured 
through the Doppler shift of the buoy signal received by the satellite, to be within 
300 m for two thirds of the reported positions. Very few fixes were reported with 
errors greater than 10 km and these can be easily removed from the record by 
careful editing. 
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Figure 1 Schematic illustration of differences between (a) first year ice and 
(b) multiyear ice (after Comiso, 1983). 



16 




SSM/I Frequencies 



4 5 6 7 8 9 10 20 30 40 100 




1 0 



0 9 



0 8 



07 



0 6 



05 



0 4 



0 3 



0 2 



Figure 2 Emissivities of calm seawater, multiyear ice and first year ice as a 
function of frequency for horizontal and vertical polarizations 
(modified from Svendsen et al M 1983). 
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Figure 3 Illustration of the four terms in the radiation model. 
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Figure 4 Multispectral analysis of brightness temperatures (from Comiso, 
1983). 
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Figure 5 SSM/I north polar grid (from NSIDC, 1990). 
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m. CASE STUDY: ANOMALOUS DIURNAL CURRENTS ON 

THE YERMAK PLATEAU 



A. BACKGROUND 

1. The Yermak Plateau 

The Yermak Plateau (Figure 6) is a hooked, finger-like extension of the 
Barents Sea continental slope, located off the northwestern coast of Svalbard. 
The plateau extends 400 km into the Arctic Ocean and covers roughly 20,000 
km^. Plateau depths range from 1000 to 600 m. The shallowest depth is 425 m 
on a ridge striking north-northwest near the plateau's central eastern edge (Perry, 
1986). The plateau slopes down rapidly to 3500 m to the north and northwest. 
The eastern slopes gradually level out to a semi-enclosed 2000 m deep basin. 

2. Regional Tides 

Tides in the Arctic Ocean are predominantly semi-diurnal. Sclnviderski 
(1986) developed a model for global tidal prediction based on a 1° by 1° grid 
with accuracies better than 10 cm. Figures 7 and 8 show the structure of the 
leading semi-diurnal (M2) and leading diurnal (Kl) partial tides, respectively, for 
the Nordic Seas. Dashed lines depict the co-range lines measured in centimeters. 
The M2 amplitude varies from zero at the aniphidromic points to 80 cm and 
ranges from 20-30 cm over the Yermak Plateau. The amplitude of the Kl tide 
varies only from zero to 10 cm, ranging from 2-5 cm over the Plateau. 

Gjevik and Straume (1989) developed a regional tidal model to simulate 
the M 2 and Ki tides for the Arctic Ocean and Nordic Seas with a 50 by 50 km 
grid size. The M 2 and Kl amplitudes are similar to Schwiderski's global model. 
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The regional model also shows clockwise Kl current rotation along the west 
Norwegian shelf up to Svalbard. These authors believe this indicates the 
presence of a coastal trapped wave modal structure. 

3. Observations 

Anomalous diurnal currents over the Yermak Plateau were first reported 
by Hunkins (1986) from FRAM III and IV current meter data. The meters drifted 
from deep (4000 m) water to shallow water over the plateau's northwestern flank. 
Hunkins resolved absolute currents into orthogonal components in the along 
slope (60°) and cross slope (150°) directions and found that strong diurnal 
currents exist along the plateau's slope. Figure 9 shows current velocities at 104 
m along with the bathymetric profile for the FRAM IV drift. During the deep 
water portion of the drift, the current record is noisy with some evidence of semi- 
diurnal tidal motion. Over the slope diurnal motion dominates the signal. Once 
on top of the plateau, the semi-diurnal signal returns. Kinetic energy spectra 
confirm the dominance of diurnal over semi-diurnal amplitudes (Figure 10). 

Tidal motion polarization differs for the two tidal periods. Figure 1 1 
shows a broad diurnal peak in the kinetic energy spectra near 150° and a 
minimum near 60°. Conversely, the semi-diurnal motions are nearly parallel to the 
slope with a maximum at 60° and a minimum at 150°. The semi-diurnal alignment 
agrees with the modeled M 2 tides for the Arctic Ocean (Kowlik, 1979 as cited by 
Hunkins, 1986). Hunkins believed the cross-shelf diurnal currents indicated 
strong topographic vorticity waves. 

Strong diurnal currents were also reported from MIZEX 84 data. 
Morison et al. (1987) noted semi-diurnal currents dominated over most of the 
current meter drift record with amplitudes of 10 cm/s. Diurnal motions, however. 
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predominated when the current meter drifted over the southwestern slope of the 
Yermak Plateau. The diurnal component of motion ranged up to 12 cm/s. 

The observations of Hunkins (1986) and Morison et al. (1987) were 
confirmed by Hoffman's (1990) analysis of ARCTEMIZ 86 drifter buoy tracks. 
He showed a strong correlation between the enhanced diurnal currents and 
bathymetry. While in deep water no significant peak in the energy spectrum at 
the semi-diurnal or diurnal period was found. However, once the buoys crossed 
the 3000 m contour into shallower water the energy spectra exhibited a diurnal 
peak which exceeded the semi-diurnal peak. Figure 12 is a buoy track that 
clearly shows diurnal clockwise looping over the plateau's northern slope. 
Hoffman also reported a periodic increase in diurnal clockwise motion that he 
associated with the fortnightly (spring-neap) beat in surface height inequality 
caused by the superposition of the two principal semi-diurnal tidal components. 

4. Analysis 

Hunkins (1986) showed that topographically trapped vorticity waves 
can exist at the diurnal frequency over the Yermak Plateau. He ignored the 
details of forcing and used a simple free wave analytic model. Chapman (1989) 
explored the actual forcing of the waves by the diurnal tide. He showed that a 
rectilinear oscillating barotropic current is significantly amplified by Hunkin's 
idealized representation of Yermak Plateau topography. 

Recent work by Prazuck (1991) indicated high kinetic energy cells are 
possible on the Plateau. His linear barotropic numerical model created such cells 
in close agreement with observed diurnal currents. He showed through an 
analytic model that the highest energy cell could be due to the backscattered 
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energy of a coastally trapped wave that is unable to propagate beyond the tip of 
the plateau. 

Prazuck concluded that the enhanced rotating diurnal currents may 
influence ice floe motion. The high energy cells would increase the grinding 
together of the floes and alter the break-up and freeze-up patterns in the region. 
Several phenomena might be related to this interaction: 

• Strong turbulent fluxes associated with the diurnal currents 
such as those found during CEAREX 89 by Padman and Dillon 
(1991). 

• Enhanced ambient noise in the region of grinding. Cousins 
(1991) correlated strong variations in ambient noise levels with 
12 hour tidal/inertial oscillations east of the Yermak Plateau. 

• Altered ice morphology due to the convergence/divergence 
of ice floes. 

5. Effects on Sea Ice Characteristics 

The interaction of the enhanced diurnal currents on the sea ice should 
be a seasonal phenomena. The high ice concentrations of the winter pack ice 
should significantly dampen floe motion, reducing the possibility of altering ice 
conditions. Lower ice concentrations in summer on the other hand should permit 
larger amplitude motions and increase the chances of remotely observable 
changes of sea ice characteristics. 

Ice floe motion effects both ice concentration and morphology. The 
enhanced currents can induce local divergence in the ice field and raise or lower 
ice concentrations on a scale observable from space. Divergence and grinding 
together of floes can increase fracturing. Subsequent refreezing would increase 
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the fraction of first year ice in a predominantly multiyear ice field without 
significantly affecting the total ice concentration. 

B. METHODOLOGY 

Data from the NSIDC SSM/I CD-ROM data set were used to test the 
hypothesis that changes in ice concentration/morphology could be detected by 
satellite remote sensing. SSM/I data were used because, although visible and IR 
measurements provide higher resolution data, cloud contamination impede most 
attempts to examine a time series extending beyond a few days. 

Sea ice concentrations were computed using the algorithm adopted by the 
NASA Sea Ice Algorithm Working Group (NSIDC, 1990). Gradient ratios were 
also calculated by modifying the ice concentration program to return gradient 
ratios instead of ice concentrations. The 25 km grid spacing used by the NSIDC 
data set was retained. Appendix A contains the software code listing. 

Figure 13 shows how the Yermak Plateau region was divided into test and 
control areas. A strip of 5 pixels (125 by 25 km), extending over the 
northwestern flank of the plateau, was selected as the test area. Control areas of 
the same size were located 50 km to the northwest and 125 km to the southeast. 
These strips will be referred to as 'Slope', 'North' and 'South', respectively. A 
perpendicular strip 75 km upstream of the general ice flow, northeast of the 
plateau tip, was used to establish the gradient ratio of ice flowing into the test 
area. 

A 168 day time series of daily ice concentrations, extending from 1 May to 15 
October 1988, was studied. One four day and several one to two day data gaps 
in the record were filled by linear interpolation of surrounding values. As seen in 
Figure 14, ice concentrations varied from 75-100% during the period for the 
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northern and test strips, and dipped to as low as 35% for the southern strip. At 
other times of the year concentrations remained near 100%. For the gradient ratio 
analysis the time series extended from 13 January to 30 April 1988. This period 
avoided brightness temperatures contaminated by melt ponds. 

C. RESULTS 

a. Ice Concentration 

The ice concentration algorithm was first checked for spatial 
consistency. Figure 15 shows concentration values for each pixel within the 
Slope strip during the period of interest. Table 3 contains a correlation matrix 
showing that the values are spatially consistent with 0.74 as the lowest 
correlation between the two endpoint pixels. 

Differencing between the North - Slope and South - Slope 
concentration values (Figure 16) showed no significant deviations between ice 
concentrations over the slope and control areas. The North - Slope mean 
difference is 1.04%, well inside the 10% accuracy of the concentration algorithm. 
The Slope - South mean difference is larger at 8.56%, but as seen from the 
concentration record this is due to the generally lower concentrations in the 
South strip later in the summer. 

The record was next examined for evidence of the fortnightly beat of 
current intensity reported by Hoffman (1986). If the currents influence ice 
concentrations, the largest signal should be present when the currents are 
strongest. An autocorrelation of Slope concentration values showed no 
significant correlations over a fortnightly period (Figure 17). Spectral densities 
were also computed for the mean Slope, North and South concentrations, as well 
as the North - Slope and South - Slope differences (Figures 18-19). In none of 
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these cases was a spectral peak located at or near the 13.66 day lunar fortnightly 
cycle (0.07 cycles per day). 
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b. Gradient Ratio 

The gradient ratio is indicative of ice age. As shown in Table 2, the 
gradient ratio should be near zero for first year ice, less than zero for multiyear ice 
and greater than zero for sea water. The Siberian Branch of the Polar Drift Stream 
advects mostly multiyear ice through this region (Polarstern Scientific Party, 
1988). If interactions between the diurnal currents and the ice are sufficiently 
strong to increase ice fracturing over the Yermak Plateau slope, the refreezing of 
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the resultant small leads should increase the total fraction of first year ice. An 
increased fraction of more saline ice will drive the gradient ratio more positive 
compared with the control areas to north, south and upstream of the flow. 

Gradient ratios for the North, Slope and South regions are shown in 
Figure 20. Mean values are -0.033, -0.031 and -0.034, respectively. All gradient 
ratios are negative as expected from the predominantly multiyear ice. Although 
the gradient ratios are at least 26% higher over the slope than over the control 
areas during days 10-42, the difference between upstream and slope values is 
negligible for the same period (Figure 21). This means that no change measurable 
by the SSM/I occurred as the ice flowed over the slope region. 

D. CONCLUSIONS 

Anomalous diurnal currents over the Yermak Plateau slopes have been 
reported and analyzed by Hunkins (1986), Morison et al. (1987), Hoffman (1990) 
and Prazuck (1991). These currents have the potential for affecting the sea ice 
cover in this region in several ways. Divergence/convergence patterns could 
alter local ice concentrations. Stresses induced in the ice could lead to increased 
grinding and fracturing of floes, changing the morphology of the ice cover. If 
such changes could be noted through remote sensing, existing data bases could 
be explored for similar phenomena in other ice-covered regions. 

Ice concentration/morphology changes are not observed by the SSM/I. Two 
explanations could account for this lack of confirmation: 

1 ) The currents have little or no effect on the ice cover. 

2) SSM/I limitations obscure features in the ice caused by the currents. 
Factors supporting the second hypothesis include poor SSM/I spatial resolution 
and contamination by melt ponds or snow cover. Passive microwave emissions 
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observed by the SSM/I are too imprecise to verify the anomalous current effects 
on the sea ice cover. 
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Figure 6 Bathymetric chart of the Yermak Plateau (from Hoffman, 1990). 
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Figure 7 Leading semi-diurnal (M 2 ) ocean tide (from Schwiderski, 1986). 
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Figure 8 Leading diurnal (Kj) ocean tide (from Schwiderski, 1986). 
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Figure 9 




Absolute current velocities along the FRAM IV track at 104 m (from 
Hunkins, 1986). 
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Figure 10 Kinetic energy spectra for the cross slope velocity component from 
FRAM IV at 303 m (from Hunkins, 1986). 
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Figure 11 Peak amplitudes of kinetic energy for the Ki and M 2 frequencies as a 
function of current axis orientation for FRAM IV at 303 m (from 
Hunkins, 1986). 
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Figure 12 Diurnal looping of ARCTEM1Z 86 drifter buoys (from Hoffman, 
1990). 



36 




Figure 13 Yerniak Plateau test and control areas. Hatched boxes indicate 
individual pixels. Depths are in meters. 
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Figure 14 Total ice concentrations, 1 May - 15 October 1988: (a) North 
control, (b) Slope and (c) South control. 
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Figure 15 Yermak slope ice concentrations, 1 May - 15 October 1988. 
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Figure 16 Percent concentration differences (1 May - 15 October 1988): (a) 
North control minus Slope, (b) South control minus Slope. 
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Figure 17 Autocorrelation of slope region ice concentrations. 
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Figure 18 Ice concentration spectra (1 May - 15 October 1988): (a) North 
Control, (b) Slope, and (c) South control. 
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Figure 19 Ice concentration spectra (1 May - 15 October 1988): (a) North - 
Slope difference and (b) Slope - South difference. 
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Figure 20 Gradient ratios for Slope and Control areas, 13 January - 30 April 
1988. 
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Figure 21 



Slope - Upstream gradient ratio difference, 1 February - 30 April 
1988. 
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IV. CASE STUDY: ICE MOTION ANALYSIS 



A. BACKGROUND 

1 . General Ice Flow 

Sea ice motion affects a variety of important arctic processes. Carsey 
and Holt (1987) note that ice movement affects latent heat advection, oceanic 
surface stress, open water production through divergence and shear, is a passive 
tracer of ocean currents, and affects icebreaker navigation and ice structure. 
Within the central ice pack, relative ice motion creates both open water, through 
which significant heat is lost, and ridges that are operational hazards. At the ice 
margins ice growth associated with divergence affects brine production which, in 
turn, influences intermediate and deep water formation. 

Ice motion is primarily due to local surface wind stress, ocean currents, 
Coriolis and inertial forces, and effects of horizontal stresses transferred through 
the ice itself (Hibler, 1986). Colony and Thorndike (1984) examined a 90 year 
record of ship, research station and buoy tracks to create a quantitative estimate 
of the mean ice motion in the arctic basin (Figure 22). Two features dominate the 
circulation pattern — the anticyclonic Beaufort Gyre in the Canadian Basin and 
the Transpolar Drift Stream (TDS) extending across the Eurasian Basin to Fram 
Strait. 

The TDS is the largest and most concentrated meridional ice flow in the 
world. It is the primary outflow region for arctic sea ice. Variations in flow 
magnitude are theorized to effect the global climate (Vinje and Finnekasa, 1986). 
As seen in Figure 23, the TDS is composed of two branches. The Polar Branch 
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advects old multiyear ice across the pole and out through Fram Strait. The 
Siberian Branch contains thick, deformed multiyear ice floes that probably 
originated on the Siberian shelves. Closer to the Yermak Plateau, generally thin, 
first year ice formed locally or in the Barents and Kara Seas is found (Polarstem 
Scientific Party, 1988). The merging of these different ice types makes Fram Strait 
a rich area for SSM/I feature tracking. 

2. Feature Tracking 

Feature tracking estimates motion of a fluid by tracking features that are 
advected in the fluid. Zhang (1991) identifies two main steps necessary for 
successful feature tracking: 

1) Tie point selection - Identification in the first image of suitable 
features to be tracked; and 

2) Tie point registration - Matching the same tie points in a subsequent 

image. 

The underlying assumption is that the tracked features are advected by 
the fluid. Motion vectors are representative of the surrounding medium only 
when this assumption is valid. Possible disruptive influences include waves in the 
fluid, when fluctuations in fluid properties are tracked rather than the Lagrangian 
flow of the fluid, atmospheric contamination which can temporarily obscure 
features, and Eulerian effects such as local heating or cooling that provide feature 
sources and sinks. 

In the past sea ice motion has been inferred from the drift records of ice- 
locked ships, manned research stations and buoy tracks. The number and spatial 
density of these in-situ observations is limited. Satellites offer the advantage of 
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repeated sampling over large areas — a potentially rich source of information on 
the spatial and temporal structure of the sea ice motion field. 

Sea ice movement can be computed from displacements of features in 
sequential images. Subjective techniques, which require manual tracking of 
specific floes, are labor intensive and subject to repeatability problems due to 
differences in operator training. Automated techniques, such as the maximum 
cross correlation method described by Ninnis et al. (1986), are much more efficient 
but must be heavily filtered to remove spurious movement vectors. 

Three electromagnetic regimes are currently available for passive space- 
based observations of ice motion. In order of decreasing resolution, they are 
visible (0.3 - 0.7 (im), infrared (3-12 (J.m) and microwave (0.3 - 1.5 cm). The major 
limitations of high resolution visible and IR imagery are their inability to penetrate 
cloud cover and the light dependence of visual wavelengths. For example, 
Emery et al. (1991) investigated Fram Strait ice motion with IR imagery but were 
limited to six images over a 10 day period in April 1986 that contained 
exceptionally cloud-free conditions. On the other hand, low resolution passive 
microwave imagery can penetrate clouds and are independent of light conditions. 

Active microwave systems offer the advantages of cloud penetration 
and light independence with much higher resolutions than passive methods. For 
example, Synthetic Aperture Radar (SAR) data from SEASAT had a resolution of 
25 m (Carsey and Zwally, 1986). SAR information from the ERS-1 satellite will 
soon be available, however its application will be limited due to the limited 
sampling in time and space, especially in the subpolar regions (Zhang, 1991). 

Many researchers believe satellite passive microwave sensing is 
unsuitable for tracking ice motion (Zhang, 1991; Cavalieri, 1991 personal 
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communication) due to its inherent coarse spatial resolution. The primary use of 
this satellite data for arctic sea ice studies has been to examine the temporal and 
spatial variability of the ice cover. Svendsen et al. (1983) reported the motion of 
a patch of high multiyear ice concentration that correlated well with buoy 
movements, but the potential for feature tracking has not been pursued in the 
published literature (Cavalieri, 1991 personal communication). 

One can understand why passive microwave feature tracking has been 
overlooked. As mentioned earlier, brightness temperatures are sensitive to 
variations in surface conditions such as wet snow or melt ponds. The spatial 
resolution is too coarse to track individual floes. Judicious application of SSM/I 
data, however, can overcome these limitations. The basic requirements are large 
scale spatial gradients in ice type (as determined from gradient ratios) to overcome 
the coarse data resolution, and small temporal variability of surface conditions to 
limit melt pond or wet snow contamination. The first requirement is achieved by 
examining areas where old and new ice merge. The second requirement is met by 
rejecting data during the summer melt season. 

These two restrictions limit the wide-scale use of inferring ice motion 
through SSM/I feature tracks. Regions where this technique shows the most 
promise, however, are areas of high interest such as Fram Strait, through which 
occurs approximately 90% of the heat exchange and 75% of the mass exchange 
between the Arctic and World Oceans (Aagaard and Greisman, 1975). 

B. METHODOLOGY 
1 . Satellite Images 

Multiyear ice fractions were calculated using the NASA gradient ratio 
technique described earlier. The same software that computed ice concentrations 
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and gradient ratios for the Yermak Plateau case study was used to create 
multiyear ice fraction images. Appendix A contains the FORTRAN code listing. 
The software merges the computed multiyear ice fractions with a land mask to 
produce a 304 by 448 pixel polar stereographic representation of the arctic 
region, extending from the pole south to 70° N. The pixel resolution is 25 km. 

Images were produced from the NSIDC Winter 88 and Winter 89 data 
sets. Images were created at 3 day intervals. If one assumes a typical 10 cm/s ice 
speed, a feature advected by the flow will move nearly 26 km over 72 hours. 
Since the grid resolution is 25 km, a 3 day sampling interval is the minimum 
needed to show expected movement of a feature within the data set. 

The imagery was displayed on a Macintosh II desktop computer using 
Imagic image processing software (Norris and Emery, 1989). Imagic was designed 
to display satellite remote sensing imagery on the Macintosh and is distributed by 
NSIDC with the CD-ROM data set. Images were color enhanced to accentuate 
multiyear ice features and looped to provide a continuous animation of feature 
movement. 

2. Feature Tracking 

The two main steps for determining ice motion from successive images 
are tie point selection and registration. These steps are achieved in this study 
through subjective imagery analysis following rules that could be incorporated 
into an objective scheme. Tie points are selected by identifying spatially coherent 
features through a visual inspection of the animated imagery. These tie points are 
registered in succeeding images by matching the same points in their new 
locations. This procedure is described in detail below: 

1) Create a series of multiyear ice fraction data sets. 
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2) Load the data sets into Imagic and generate ice fraction images. 
Decide on a color enhancement scheme that best shows multiyear ice features 
and apply this enhancement to all images. 

3) Using Imagic, observe the advection of the multiyear ice features. 

4) Select features that appear to be advected with the general ice flow 
according to the following criteria: 

• Size - Features must be small, no more than 100 km across (4 
pixels. Accurate tiepoint registration within larger features is too 
difficult. 

• Background Discrimination - Features must stand out from the 
background by at least 10 percentage points for most of their 
lifespan. This increases one’s confidence that what is being 
observed is actually a feature advected in the flow rather than 
random variability of multiyear ice fractions. 

• Temporal Coherence - Features must be present in a least 6 
consecutive images (18 days). This further reduces the possibility 
of spurious tracks due to temporary brightness temperature 
fluctuations in the ice cover. 

5) Record the pixel position of an appropriate tie point within the 
feature during its lifespan using the following guidelines: 

• If one pixel value is 3 percentage points or more higher than all 
other values within the feature, select it. 
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• If 2 or more pixel values are within 3 percentage points, select 
the pixel which best corresponds to the previous relative position 
of the tie point within the feature. 

♦ If the relative maximum within the feature is ambiguous, 
choose the pixel which most conserves ice motion. 

6) Convert the pixel location within the polar stereographic grid to 
geodetic latitude and longitude using software supplied by NSIDC. A version 
modified for the Macintosh is listed in Appendix A. 

7) Determine ice velocities by measuring the change of position with 
time. FORTRAN code that converts geodetic positions to u and v components of 
velocity is also contained in Appendix A. 

An example of the feature tracking technique is shown in Figures 24 
through 26. This particular feature was tracked from 1 February through 22 April 
1988 as it moved southward down the east coast of Greenland. Figure 24 
contains six consecutive images of multiyear ice fraction contour plots from 3-17 
March. Each plot covers the same 500 by 250 km region centered at 76.1 °N, 
13.5°W (note the Greenland coast in the lower left corner of the plots). An area 
of relatively high multiyear ice fractions can be seen migrating from the upper left 
to the lower right comer of each plot. This is the trackable feature. Figure 25 
shows pixel values for this same region expressed as multiyear ice fractions. The 
circled number is the feature tie point. Pixel locations were converted to geodetic 
coordinates and the u and v velocity components are shown in Figure 26. 

3. Feature Selection 

The number of trackable features depends on the region of interest. For 
example, the generally homogeneous conditions in the high arctic (above 85°N) 
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create a poor environment for feature tracking. No suitable features were 
observed in the 1988 data set. Only two features in the high arctic were 
observed during the 1989 winter season. On the other hand, Fram Strait is a 
potentially rich area due to the merging of the Siberian and Polar branches of the 
TDS. 

To assess the validity of this feature tracking technique, the trajectories 
of multiyear ice features were compared with true ice motion using the arctic data 
buoy program as ground truth for the ice velocities. The criterion to colocate 
trackable features with a drifting buoy was challenging. Suitable features had to 
be identified and their positions compared against buoy tracks as reported in 
Colony and Rigor (1990 and 1991). Three feature-buoy matchups were selected 
from the winter 88 data set, each representing a different region of the arctic. 
Named after the closest geographic reference point, the features are Severnaya 
Zemyla (SZ), Franz Josef Land (FJL), and East Siberian Sea (ESS). A feature from 
the Winter 89 data set (CA) was chosen to represent motion in the central arctic 
basin. Although no satisfactory matchup between a trackable feature and a 
drifting buoy was found in the Fram Strait region, a buoy (#7001) that passed 
down the Greenland coast two weeks ahead of a feature (FRAM) was used for a 
relative comparison of the ice trajectory. 

The tracked features and ground truth buoys were not coincident in 
space; the average distance between the selected features and buoys varied from 
76 to 557 km. The closest matchup (76 km) occurred between FJL and buoy 
#1790. The ESS feature was 420 km from buoy #1896. SZ was 320 km from the 
buoy #7048. CA averaged 497 km from the closest buoy (#3288). Table 4 
summarizes the characteristics of the feature-buoy tracks. 
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TABLE 4: MYI FEATURE - BUOY DATA 



Feature 


Buoy 


Separation 

(km) 


Start Track 
Date 


Position 


Stop Track 
Date 

(Total Days) 


Position 


ESS 


7048 


420 


1 APR 88 


76.2N, 
160. IE 


25 APR 88 
(24) 


77.3N, 

159.0E 


FJL 


1790 


76 


1 FEB 88 


80.2N, 

41.3E 


14 MAR 88 
(42) 


79.5N, 

41.9E 


SZ 


1897 


320 


1 FEB 88 


81. 7N, 
105.9E 


28 APR 88 
(85) 


82.0N, 

100.6E 


CA 


3288 


491 


1 JAN 89 


84.5N, 

143.5E 


31 MAR 89 
(90) 


86.6N, 

125.2E 


FRAM 


N/A 


N/A 


1 FEB 88 


79.9N, 

4.6W 


22 APR 88 
(80) 


68.1N, 

23.9E 



Gcolocaiion Key: SZ - Severnaya Zemlya Note: Features are named 

FJL - Franz Josef Land according to nearby 

ESS - East Siberia Sea geographic reference 

CA - Central Arctic points. 

FRAM -Fram Strait 



4. Errors 

Before ice motion vectors between the buoy - feature pairs are 
compared, the maximum errors inherent to this technique are described. Three 
types of positional errors are possible: imperfect tiepoint registration, SSM/I data 
geo-location, and ARGOS buoy positioning. In addition, ice velocity gradients 
between the feature and buoy locations need to be considered. These errors, 
expressed as an ice speed error over a 3-day period, are examined in detail below. 

Imperfect tiepoint registration within the data grid is the largest source 
of error in this feature tracking technique. For small, well-defined features the 
maximum registration error is estimated to be one pixel for each observation. The 
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actual error should be uniformly distributed within the pixel. The maximum error 
is equivalent to a 50 km or 19.3 cm/s velocity error over a 3 day period. The geo- 
location error, as reported by NSIDC (1990), is 8 km for each scan. Since each 
pixel represents a daily averaged mean of multiple scans, the geo-location error 
should be somewhat smaller. For an error of 5 km, the magnitude of ice motion 
could vary by 3.8 cm/s over 3 days. Thus the maximum estimated position error is 
30 km. If this error is assumed to be uniformly distributed, then the standard 
deviation of the error is 19.0 km. 

The estimated velocity error due to the maximum position error is 23.1 
cm/s over 3 days. This is over twice the typical ice velocity of 10 cm/s. If instead 
the standard deviation of the position error is considered, the velocity error 
reduces to 14.7 cm/s. These velocity errors decrease if the motion is averaged 
over longer time periods. For example, a 12 day average would result in a 
maximum error of 5.8 cm/s. 

The ARGOS system currently used to track buoy locations can 
determine position accuracy within a few hundred meters (Thorndike, 1986). 
Assuming a maximum error of 500 m for each ARGOS reported position, the 
resulting velocity error over 3 days would be 0.4 cm/s. 

When comparing SSM/I-derived ice trajectories with buoy locations, the 
velocity gradient between the two locations must be considered. These velocity 
gradients are the result of the ice pack response to external forces such as wind 
and ocean currents. Thorndike (1986) compared the velocity components of ice 
camps and buoys separated by distances ranging from 1 to 1000 km. He 
concluded that on scales of this magnitude ice behaves as a homogeneous, 
isotropic velocity field. Thorndike examined the correlations of buoy and ice 
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camp velocities as a function of distance. For a 100 km camp-buoy separation, he 
found the correlation to be 0.9. At a 400 km separation, the correlation drops to 
approximately 0.6. 

Using these correlations, the variance of the differences in ice velocity 
between two separate points within the ice pack can be estimated. The variance 
of the velocity difference between the points (cr 2 e rr) is given by 

a 2 err = E[V(x+s)-V(x)] 2 (8) 

where E is the expected value operator and V(x) and V(x+s) are the ice velocities 
separated by some distance s. Assuming the variance of ice velocity is the same 
at both points, ( 8) reduces to 



0 2 err = 2a 2 i(l-p) (9) 

where o 2 i is the variance of the ice velocity and p is the correlation coefficient. 
Given a standard deviation of ice velocity of 10 cni/s (Thorndike, 1986) and 
assuming the u and v components of motion are independent of each other (the 
total variance then equals the sum of the variance of each component), o 2 err is 40 

cm 2 /s 2 and 10 cm 2 /s 2 for correlations of 0.6 and 0.9, respectively. The resulting 
velocity errors are 6.3 cm/s for a 400 km separation and 3.2 cm/s for a 100 km 
separation. These velocity errors equate to 16 km and 8 km positional errors, 
respectively. 

Table 5 summarizes the maximum errors expected for each error source. 
The first three position errors are per observation. The corresponding velocity 
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error results from combining two position errors for each velocity calculation. 
Velocity errors are calculated over 3 and 12 day day intervals to show how the 
maximum estimated error decreases for longer mean flow calculations. 



TABLE 5: SUMMARY OF MAXIMUM ERRORS 



Error Source 


Positional Error 
(km) 


3 - Day 
Velocity Error 
(cm/s) 


12- Day 
Velocity Error 
(cm/s) 


Tie Point 
Registration 


25 


19.3 


4.8 


SSM/I 

Geolocation 


5 


3.8 


0.9 


Buoy 

Position 


0.5 


0.4 


0.1 


Velocity 
Gradient 
(100 Km) 


8 (3 Day) 


3.2 


3.2 


Velocity 
Gradient 
(400 Km) 


16 (3 Day) 


6.3 


6.3 



This error analysis shows that 3 day ice motion estimates are unreliable 
indicators of true ice motion because the velocity error can exceed the typical 
true ice velocities of 10 cm/s. Longer term averages, however, reduce the impact 
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of the SSM/I position errors and have the potential to represent accurately 
general ice flow. 

An expected correlation coefficient can be calculated from the maximum 
expected errors. This correlation can be used to validate the hypothesis that the 
feature tracks are representative of true ice motion. If the observed correlations 
are less than the expected value, then either the feature tracks do not represent 
true ice motion or the errors were underestimated. If the observed correlations are 
larger than the expected value, then the actual errors are smaller than expected. 

The correlation coefficient between two quantities, x and y, is given by 



Exy 

[Exx Eyy]^/“ 



( 10 ) 



where both x and y are assumed to have zero means. Let x and y be defined as 
the ice motion at the feature and buoy, respectively, as observed from the data: 

x = v+ef+ed (11) 

y = v + eb (12) 



where v is the true ice velocity, ef is the error in measuring the velocity of the 
feature, ed is the error due to the velocity difference between the feature and the 
buoy and eb is the error in measuring the velocity of the buoy. The variances due 
to these errors are defined as 



G“f = Eef2 + Eed- 
o-b = Eeb 2 



(13) 

(14) 
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(15) 



o 2 i = Ev 2 



where c 2 f is the variance caused by the feature position and velocity difference 
errors, a 2 b is the variance due to the buoy position errors and a 2 i is the variance 
of true ice velocity. The expected values in equation (10) become 



E X y = o 2 i 
Exx = a 2 i + a 2 f 
Eyy = 0 2 i + C 2 b 



(16) 

(17) 

(18) 



Thus the expected correlation coefficient can be defined as the variance 
of true ice motion divided by the square root of the variances of buoy motion and 
feature motion multiplied together: 



a 2 ; 



P = 



[(a 2 ; + a 2 b )(G 2 i + o 2 f )]l/ 2 



(19) 



As before, a 2 j for each component of motion is 50 cm 2 /s 2 . Using the 3 day 
sampling interval from Table 5, o 2 b is negligible at 0.16 cm 2 /s 2 and o 2 f is the sum 

of the variances due to positional error (266.8 cm 2 /s 2 ) and the velocity gradient 
between the buoy and the feature (39.7 cm 2 /s 2 for 400 km and 10.2 cm 2 /s 2 for 
100 km separations). Using these values in equation (19), the expected 
correlation for a feature and buoy 100 km apart is 0.391. For a separation of 400 
km, the correlation decreases to 0.374. 
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These numbers indicate that if the correlation between 3 day mean 
velocities determined from buoys and multiyear ice features exceeds 
approximately 0.4, the feature is moving with the ice and the errors associated 
with the tracking technique were less than expected. If the same calculations are 
performed for 12 day mean ice velocities, correlations rise to 0.866 and 0.783 for 
100 and 400 km separations, respectively. 

C. RESULTS 

1. Ice Trajectories 

Trajectories of the four feature - buoy pairs from the different arctic 
regions, each plotted at 3 day intervals, are shown in Figures 27 through 31. 
Because the feature tracks are limited to discrete 25 km jumps in position, periods 
of slow, looping or retrograde motion confuse the trajectories. In these cases, the 
positions have been slightly offset to show movement more clearly. 

The correspondence between 3 day positions of the buoys and features 
appears somewhat erratic. This is probably due to the resolution limitations 
inherent in the SSM/I data grid. Taken over longer periods, however, the general 
flow corresponds quite well to the buoy positions. This is best shown in the Fram 
Strait trajectories (Figure 27). Both the buoy and FRAM tracks are roughly 
parallel as they traveled southward down the east Greenland coast, the buoy 
preceding FRAM by about two weeks. Of special note is the sudden surge 
eastward at 78°N for the feature and at 69°N for the buoy. Although separated 
by 9° of latitude, this surge occurred at the same time (20-23 Feb). This 
corresponds to a period of strong geostrophic wind flow to the northeast due to a 
984 mb low centered over the northern coast of Greenland (Colony and Rigor, 
1990). 
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2. Velocity Comparisons 

For each feature - buoy pair, the u and v components of motion over 
each 3 day period were calculated. Buoy velocities were computed directly from 
their reported positions. With the exception of ESS, all multiyear feature 
velocities were smoothed over two positions to reduce velocity spikes resulting 
from the gridded data set. The ESS feature, tracked for 24 days, has only one 
local velocity maximum in each direction so the smoothing was not applied. 

Correlation coefficients between the buoys and both smoothed and 
unsmoothed feature velocities are presented in Table 6. Applying a t distribution 
test to the smoothed velocities, all but one correlation (Franz Josef Land u 
component) exceed the 95% confidence level that the feature and buoy 
velocities were correlated. The generally high correlations are particularly 
surprising since the maximum error analysis indicates that correlations only as 
high as 0.4 should be expected for 3 day mean ice velocities. These higher than 
expected correlations imply that the estimated errors used in the error analysis 
were too large. 

Specific velocity comparisons are plotted in Figures 32 through 35 and 
are discussed below. 

a. Franz Josef Land - Buoy 1790 (Figure 32) 

This is the closest spatial match-up of any feature - buoy pair 
investigated. The separation distance averaged 76 km over a 42 day period. The 
buoy trajectory shows complex motion with two separate periods of retrograde 
east - west movements. The SSM/I grid had trouble resolving these complicated 
movements, providing the lowest correlation of any velocity component: -0.069 
for the unsmoothed and 0.322 for the smoothed u component. When velocities 
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were averaged over 12 days, the unsmoothed u component correlation improved 
to 0.666. Three day v correlations were 0.458 and 0.678 for unsmoothed and 
smoothed velocities, respectively. The higher correlations in the v direction are 
probably due to the slow, southward trend of the drift. 



TABLE 6: MYI FEATURE - BUOY CORRELATIONS 



Feature 


Buoy 


Correlations (3 day velocities) 


u 


V 


smooth u 


smooth v 


ESS 


1896 


.791 


.724 


.845 


.490 


FJL 


1790 


-.069 


.458 


.322 


.674 


SZ 


7048 


.652 


.537 


.685 


.455 


CA 


3288 


.051 


.501 


.311 


.582 



b. East Siberian Sea - Buoy 1896 ( Figure 33) 

This pair was tracked over the shortest interval (24 days). Although 
separated by 320 km, it also had the highest correlations. The correlations for the 
unsmoothed u and v components were 0.791 and 0.724, respectively. The 
trajectory characteristics responsible for these high correlations were a steady 
northwestern drift interrupted by a strong southwesterly surge (16-19 April) large 
enough to be resolved within the SSM/I grid. This surge corresponded to a 
period of strong geostrophic wind flow to the southwest from a low in the 
Chukchi Sea region (Colony and Rigor, 1990). 
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c. Severnaya Zemyla - Buoy 7048 (Figure 34) 

The feature and buoy were tracked for 75 days. The average 



separation distance was 420 km. Both trajectories showed a slow easterly drift 
culminating in a clockwise loop. Correlations for the smoothed u and v 
components were 0.685 and 0.455, respectively. The higher u correlation is due 
to the predominantly eastward movement of the feature and buoy. 
d. Central Arctic - Buoy 3288 ( Figure 35) 

This pair was taken from the 1989 data set to represent ice motion in 
the central arctic region. They were tracked for 90 days at an average separation 
of 491 km. The smoothed u correlation was low at 0.31 1. The v component was 
much higher at 0.582. The low u correlation can be explained by the steady 
northerly movement of the feature. East-west velocities were generally small and 
not well resolved in the SSM/I data grid. Correlations for the 12 day mean 
velocities jumped to 0.783 and 0.776 for u and v components, respectively. The 
longer sampling period allowed time for the slow motion of the ice field in this 
region to be resolved within the SSM/I grid. 

3. Position Error 

The maximum estimated feature position error can be compared against 
the position errors calculated from the data. Solving equation (19) for o-f: 



negligible compared with G~ice . The average variance of each buoy velocity 




( 20 ) 



i is assumed to be the same as the variance of buoy velocities and G-^b is 
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component is 30.08 cm 2/s 2 and the average correlation (p) between the 
smoothed feature and buoy velocities is 0.545. 

Using equation (20) and the average variances and correlations 
observed from the feature track data, c 2 f is 71.19 cm 2 /s 2 . The ice velocity 
difference between the buoy and feature must be taken into account because c 2 f 
is the sum of the variances due to position error and the velocity difference. 
Substituting the average standard deviation of the buoy velocities into equation 
(9) and again assuming that u and v components of motion are independent of 
each other, the variance due to the ice velocity gradient is 12.32 cm 2 /s 2 and 3.08 
cm 2 /s 2 for 400 km and 100 km feature-buoy separations, respectively. Taking 
into account the average feature-buoy separation was 326 km, nearly 17% of the 
observed variance can be explained by the ice velocity gradient between the 
feature and the buoy. The remaining variance of approximately 58 cm 2 /s 2 
equates to a 10.7 cm/s total velocity error (u and v velocity variances added 
together). 

This velocity error is less than half of the maximum expected value of 
23.1 cm/s and 72% of the standard deviation of the expected error. A 10.7 cm/s 
velocity error over a three day period implies close to a 15 km position error in 
each of the two observations used to calculate the velocity. This is less than both 
the maximum expected postion error (30 km) and the standard deviation of this 
error (19 km). Possible explanations for the difference between expected and 
actual values include: 1) the tiepoint registration and geo-location errors are 
smaller than expected because the gridded brightness temperatures in the dataset 
are obtained by averaging information gathered from several satellite passes each 
day; 2) the actual ice velocity gradient between the features and buoys is smaller 
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than Thorndike's (1986) estimates; and 3) the standard deviation of the ice 
velocities is smaller than the assumed value of 10 cm/s. The latter explanation 
probably had the greatest influence because standard deviations of the buoy 
velocities only ranged from 3 to 9 cm/s. These smaller than expected standard 
deviations could result from the the use of three day ice velocity averages which 
reduce the effect of short term velocity surges. 

I). DISCUSSION 

There is good reason to doubt the practicality of inferring sea ice motion by 
tracking multiyear ice concentrations from space-based passive microwave 
radiometers. Data resolution is poor. At 37 and 19.35 GHz (the frequencies 
currently used to calculate ice fractions), resolutions are 33 and 56 km, 
respectively. The 25 km grid resolution of the SSM/I dataset is only achieved by 
averaging several observations over the same spot each day. Not only does the 
low resolution increase the position error, it also greatly reduces the possibility of 
tracking individual floes. The radiometer must instead observe broad regions of 
anomalous radiances to create trackable features. The precision of the multiyear 
ice fraction algorithm is only estimated to be in the range of 13 - 25 %. Moreover, 
a variation of 0.04 cm in cloud liquid water content will contribute as much as 
30% to the uncertainty of computed multiyear ice fractions (Cavalieri et al., 
1984). 

In spite of these handicaps, ice drift inferred through tracking multiyear ice 
features from SSM/I radiometer data corresponded quite well with ice motion as 
measured by in-situ buoys. Ice feature trajectories were similar to nearby buoy 
drift patterns. In two cases (ESS and FRAM), when high atmospheric pressure 
gradients caused anomalously strong wind stress on the ice pack, both the feature 
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and buoy exhibited similar trajectories in accordance with the forcing. 
Comparisons between the orthogonal velocity components of the feature and 
buoy motion averaged over 3 days also showed surprisingly high correlations. 
Instances of high correlations in one component direction but low correlations in 
the other can be explained by slow or oscillatory movement in the direction of 
the velocity component exhibiting the low correlation and are a function of the 
low spatial resolution of the grid. 

The good results achieved above were the result of being extremely selective 
during the tiepoint selection process. First, all tracked features were of high 
multiyear ice fractions surrounded by a background of relatively low fractions. 
Multiyear ice is thicker than first year ice and therefore less subject to 
deformational and destructive processes that might obscure the characteristic 
radiometric signal. Second, the features were small — no more than 100 km across. 
This simplified the task of registering the same tiepoint in successive images. The 
variation of multiyear ice fraction percentages within a larger feature could 
preclude accurate tiepoint registration. Limiting the feature spatial extent also 
reduced the possibility of confusing an atmospheric disturbance with an ice 
feature. Third, imposing a temporal coherence requirement of 18 days ensured 
that the tracked object was not a result of random brightness temperature 
fluctuations in the ice cover. Finally, limiting the record to the winter season 
reduced the possibility of surface melt conditions contaminating radiances from 
the sea ice. Another benefit of the winter season is that moisture-laden clouds 
with their contaminating effects are generally restricted to the open ocean and 
the marginal ice zone areas. 
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The number of trackable features was limited. The abovementioned 
limitations restricted the number of suitable features to only a dozen or so during 
the winter season. Relaxing these restrictions should produce more trackable 
features (and possibly introduce spurious features) and is an area suitable for 
further research. Development of an objective scheme to identify and track" 
features should also increase the number of available datapoints. 

Despite the above restrictions, this limited ability to observe and track areas 
of anomalous multiyear ice fractions can be useful for arctic studies. Ice fluxes 
from the central arctic have been shown to be the largest exporter of fresh water 
from the Arctic Ocean and the largest freshwater input to the Greenland Sea 
(Aagaard and Carmack, 1989). Available databases, extending back to 1978, can 
be used to assess the intra and interannual variability of ice transport. 
Researchers can identify previous, and possibly forecast future, pulses of ice 
types to examine their influence on regional heat and salinity budgets. Ice camps 
are generally located on relatively smooth, thick floes of multiyear ice for 
logistical and safety reasons. Potential sites are currently scouted by aircraft. 
This space-based technique offers the possibility of remotely identifying sites well 
in advance of deployment. As satellite instrumentation improves and emanations 
from higher microwave frequencies are utilized, surface resolution will also 
improve resulting in increased reliability of short term ice velocity vectors. 

E. CONCLUSIONS 

When images of arctic multiyear ice fractions calculated from winter season 
SSM/I data are examined in sequence, distinct features are observed to move in 
the general direction of climatological ice flow. These features exist throughout 
the arctic basin but are mainly located in regions where first year and multiyear 
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ice merge together. Some features are short-lived, others can be tracked 
throughout the data record. 

Four of these multiyear ice features from the winters of 1988 and 1989 were 
compared against nearby arctic buoy trajectories. Three day average ice velocity 
vectors were found to correlate well with the buoy velocities; the average 
smoothed velocity correlation was 0.545. The maximum estimated velocity errors, 
mainly due to the low data resolution, are double a typical ice velocity of 10 cm/s 
making short-term ice motion calculations unreliable. Actual velocity errors 
derived from the feature-buoy correlations are half the maximum estimated error 
and less than three quarters of the standard deviation of this error. This difference 
is probably due to the smaller than estimated variances of ice motion because the 
ice motion is averaged over a three day period, smoothing out short-term ice 
velocity surges. 

The higher than expected correlations between multiyear ice feature 
velocities and arctic buoys are due to the stringent requirements established for 
feature selection. Relaxing these requirements could substantially increase the 
number of trackable features from the dozen or so noted in the winter of 1988 but 
also decrease the confidence that tracked features accurately represent ice 
motion. A similar increase in trackable features could probably be attained by 
objectively scanning the dataset for subtle features that were not recognized 
through a subjective analysis of the data. 
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igure 22 Field of mean ice motion (from Colony and Thorndike, 1984). 
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Figure 23 Average sea ice drift showing the two components of the 
Transpolar Drift Stream ( from Polarstem Scientific Party, 1988). 
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Figure 24a Contour plots showing advection of multiyear ice feature along the 
east Greenland coast. Images from 2 and 5 March, 1988. Plots are 
centered at 76.1 °N, 13.5°W. The feature moves from the upper left 
to the lower right corner. 
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Figure 24b Same as Figure 24a except images are from 8 and 1 1 March, 1988. 
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Figure 24c Same as Figure 24a except images are from 14 and 17 March, 1988. 
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5 March 1988 



Figure 25a Pixel values for multiyear ice contour plot shown in 
value is tiepoint registration location. 



Figure 24a. 
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Figure 25b 



Pixel values for multiyear ice contour plot shown in Figure 24b. 
value is tiepoint registration location. 
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Figure 25c 



Pixel values for multiyear ice contour plot shown in Figure 24c. 
Circled value is tiepoint registration location. 
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Figure 26 Ice velocities for Fram Strait feature (1 February - 22 April, 1988). 
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Figure 27 Ice trajectories: Fram Strait feature and Buoy 7001 (1 February - 22 
April, 1988). Circled points are initial positions. 
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Figure 28 Ice trajectories: Franz Josef Land feature and Buoy 1790 (1 

February. - 14 March, 1988). Circled points are initial positions. 
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Figure 29 Ice trajectories: East Siberia Sea feature and Buoy 1896 (1 - 25 
April, 1988). Circled points are initial positions. 



80 




Figure 30 ice trajectories: Severnaya Zemyla feature and Buoy 7048 (1 

February - 28 April, 1988). Circled points are initial positions. 
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Figure 31 Ice trajectories: Central Arctic feature and Buoy 3288 (1 January - 
31 March, 1988). Circled points are initial positions. 
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Figure 32 Franz Josef Land and Buoy 1790 velocity comparison. Buoy is 
solid line, feature is dashed line. 




Figure 33 East Siberian Sea and Buoy 1896 velocity comparison. Buoy is 
solid line, feature is dashed line. 
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Figure 34 Severnaya Zemyla and Buoy 7048 velocity comparison. Buoy is 
solid line, feature is dashed line. 
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Figure 35 Central Arctic and Buoy 3288 velocity comparison. Buoy is solid 
line, feature is dashed line. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



Arctic sea ice conditions are an important component of the global climate 
system. The sea ice cover affects the Earth's radiative balance by influencing the 
heat flux between the ocean and atmosphere and heat transport between the 
lower and upper latitudes. Passive microwave remote sensing of sea ice offers 
continuous, all-weather observations of these remote areas independent of solar 
illumination. This is particularly important in the arctic where much of the region 
is shrouded by clouds or covered in darkness. 

Accurate ice observations using passive microwave remote sensing 
techniques require knowledge of the radiative characteristics of sea ice. The 
radiative properties of sea ice are complicated due to its complex nature which is 
determined by formation environment, age and season. Other complicating 
factors include contamination from moisture-laden clouds and the presence of 
snow cover, melt ponds and ridging. In spite of these difficulties, emissivity 
differences between sea ice and water at various microwave frequencies can be 
exploited to estimate ice characteristics. Current algorithms use polarization and 
gradient ratio differences to estimate ice concentrations and type. 

Historical SSM/I data are readily available on CD-ROM media distributed by 
NSIDC. Two applications of this dataset to arctic sea ice studies were examined. 
Data from 1 February to 15 October 1988 were used to associate anomalous 
diurnal tidal currents over the Yermak Plateau with ice morphology. Data front 1 
February to 28 April 1988 and 1 January to 31 March 1989 were used to explore 
the possibility of determining sea ice motion by tracking advection of multiyear 
ice features. Software supplied by NSIDC were modified to compute and display 
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the requisite ice concentration, gradient ratio and multiyear ice fraction 
calculations and imagery on a Macintosh II desktop computer. 

The hypothesized ice concentration/morphology changes induced by the 
anomalously strong diurnal tidal currents observed over the Yermak Plateau 
could not be verified by the SSM/I data. Either the currents do not affect the ice 
cover or SSM/I limitations obscure the ice signatures of the currents. The low 
spatial resolution and susceptibility to melt and snow contamination of passive 
microwave emissions prevent a definitive assessment of tidal current effects on ice 
cover in this region. 

Distinct multiyear ice features are observed to move in the general direction 
of the mean ice flow. Feature movement correlates well with ice motion as 
measured by arctic data buoys. The subjective nature of the feature tracking 
technique and its reliability on 'features of opportunity' prevent wide scale 
application throughout the arctic region. Positional errors resulting from low data 
resolution could make short term ice motion calculations unreliable. Long term 
averages, however, are reliable and can be used to augment our knowledge of 
general ice flow. For example, motion vectors can be used to estimate the 
interannual variability of ice transport through the Fram Strait and aid our 
understanding of heat and salinity budgets in the Greenland Sea. In spite of the 
limitations, tracking multiyear ice features can be a useful method of discerning 
continuous, all-weather ice motion from satellite data. 

Further attempts to verify anomalous diurnal currents over the Yermak 
Plateau using SSM/I data are not recommended. The data resolution is too low to 
observe small scale changes of sea ice concentration or morphology. Feature 
tracking, however, does show some promise. Conversion of the subjective 
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technique described here to an objective scheme for computer analysis would 
greatly reduce the time and effort required for tie point selection and registration, 
and probably increase the number of tie points suitable for tracking. Exploration 
of the higher spatial resolutions afforded by 85.5 GHz data would greatly 
increase the reliability of ice motion vectors. 
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APPENDIX A 



This appendix contains FORTRAN source code listings of software used to 
process the dataset contained on the NSIDC CD-ROM disks. 



Program 


Paee 


Purpose 


ssmiice 


89 


Computes ice concentrations from 
SSM/I brightness temperature values on 
the NSIDC CD-ROM dataset. 


locate 


95 


Transforms I, J coordinates of the NSIDC 
SSM/I polar grids to latitude and 
longitudes coordinates and vice versa. 


u,v calc 


102 


Calculates u and v components of 
motion from a list of user specified 
geodetic latitude and longitude 
positions. 


Subroutine 


Paee 


Purpose 


usetie 


105 


Computes coefficients for the ice 
algorithm in ssmiice using NASA- 
approved tiepoints. 


ice25 


108 


Computes ice concentrations from 



NSIDC CD-ROM data at 25 km grid 
resolution. 



88 



program ssmnce 

ccccccccccccccccccccecccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c SSMIICE - This program computes ice concentrations from SSMI TB 
c grids residing on CD-ROM. 
c 

c Inputs 
c 

c USER DEFINED 



c 



c 

c drive : Name of CD-ROM disk 
c 

c indate : Date in the format YYMMDD (eg., 870710) 
c From this variable the CD-ROM is searched for the 

c appropriate filename. The filename for the ice 

c concentration algorithm is derived from this 

c variable, 

c 

c ihem : Defines the appropriate hemisphere to apply this 
c algorithm (l=North, 2=South) 

c 

c icepar : Defines the type of ice concentrations to derive: 
c l=total, 2=multi-year, 3=first-year. 

c 

c INTERNAL 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



infile : The path, name, and descriptor string describing 
the data file from which the SSMI Tb data will be 
extracted from the CD-ROM. User must specify appropriate 
filepath. 

buff : Input buffer containing the SSMI Tb data read from 
the CD-ROM. Two lines of SSMI TB grid data will 
read with a single 10 request. 

mask : Landmask used in computing the 25 km 
ice concentration grid. 

idiml : Number of parameters in an SSM/I Tb grid cell 

idim2 : Number of cells in the SSM/I Tb grid line 

idim3 : Number of grid lines in an SSM/1 Tb grid 

grdtyp : Three character descriptor designating the grid type 
to access from the CD-ROM (N3B or S3B) 
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icetyp : Three character descriptor designating the extension 
for the resultant ice parameter image (N3C or S3C). 



c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c 

c outfil : The name and descriptor string which reference 
c the data file which will be used to store the 

c ice concentration grid. This string is computed 

c internally from the variables, indate and grdtyp. The 

c file will be created in the default directory, 

c Program will abort if file currently exists! 

c 

c The file dimensions will be 304 1-byte cell values 

c with 448 grid lines for the grid type 'N3C' and 

c 316 1 -byte cell values with 332 grid lines for the 

c grid type 'S3C'. 

c 

c ice : Output buffer containing the computed ice 
c concentrations for the 25 km total ice concentration 

c grid. Totice is 3 1 6 bytes in length which corresponds 

c to the length of an S3C raster line, 

c 

c Vince Troisi 
c 5/3/1989 
c 

c 1/4/1989 - Modified for SSM/I CD-ROMs 

c 7/1/91 - Modified for Apple Macintosh and 25 km grid by W. A. Wright 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
integer idiml, idim2, idim3 
c 

c There are five channels in the N3B or S3B grid cells 
c 

parameter (idiml = 5) 
parameter (maxlen = 316) 

integer iunit, ounit, munit, olen,z 
integer ymd(6), ier, recln,ok 

character * 3 grdtyp, icetyp 
character * 9 drive 
character * 2 year, month, day 
character * 6 indate(95), ifile 
character *3 mon(12) 
character *80 outfil, infile, mfile 



c : Array of tiepoint coefficients. Returned from USETIE 

and used in GETICE. Refer to NODS Handbook Appendix 
CC pp 23-24. 

Outputs 

INTERNAL 
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integer 



* 2 buff(maxlen*idiml),abyte,bbyte 
integer * 1 ice(maxlen) , mask(maxlen),atemp(2),btemp(2) 
real c(12) 
equivalence (bbyte,btemp) 
equivalence (atemp,abyte) 
data mon 

+ / , JAN , ,'FEB', , MARVAPR , , , MAY , , , JUN , , , JUL , ) ’AUG , , , SEP , ) 
+ ’OCT','NOV','DEC'/ 



C 



c Query user for search criteria 

c Query for name of CD-Rom drive 
print *, 'Name of CD-Rom?' 
read (*, ’(a)') drive 

c Query user for file dates 

5 print *, 'Number of files to create? (95 max)' 
print *, 'Note: Read lasts approx. 10 min. per file' 
read *, n 

c print *, 'Enter file dates (yymmdd).’ 

open (unit=29, file='Archs HardDisk:THESIS:daylist', 
+ status='old') 
do 10 i=l,n 
read (29, '(a)’) indate(i) 

10 continue 

print *, 'These are the selected dates:’ 
do 20 i=l,n 
print *, indate(i) 

20 continue 

print *, 'OK? (l-yes, 2-no, 3-stop) <ret>' 
read *, ok 

if (ok.eq.2) go to 5 
if (ok.eq.3) go to 6000 

c Query user for appropriate polar region 

print *,' ' 

print *, 'Select one of the polar regions:’ 
print 1. North' 
print *,' 2. South' 
print *, 'Enter 1 or 2:' 
read *, ihem 
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c Query user for ice type to compute 



print *,' ' 

print *, 'Select type of ice characteristic to derive:' 
print 1. Total Ice Concentration' 
print *,' 2. Multi-Year Ice Concentration' 
print *,'Enter 1 or 2:' 
read *, icepar 

c Combine responses to form filenames and paths and to define 

c 10 constants 

if (ihem.eq.l) then 
grdtyp = 'N3B' 
icetyp = 'N3C' 
lrec = 3040 
idim2 - 304 
idim3 = 448 
npar = idim l*idim2 
else 

grdtyp = 'S3B' 
icetyp = 'S3C' 
lrec =3160 
idim2 = 316 
idim3 = 332 
npar = idim 1 *idim2 
endif 



do 5000 z=l,n 
ifile = indate(z) 
read (ifile(3:4),'(i2)')mm 

infile=drive//':7/grdtyp//': 19'//ifile( 1 :2)//':'//mon(mm)// 
+ ':7/ifile( 1 :6)//'. '//grdtyp//'; 1 ' 

print *, 'infile = ', infile 

c Open brightness temperature file on CD 



iunit = 21 

open (unit=iunit,file=infile,status='old',access='sequentiar, 
+ form='unformatted',recl=lrec,iostat=ier) 
if (ier.ne.0) then 

print *,infile, ' not accessible! <ret>' 
pause 
stop 
endif 
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c Construct output file name and open the file. 

outfil = ’Archs HardDisk:THESIS:Conc Data:’//ifile(l:6)// 
+ ’.'//icetyp 
print *, 'output file = outfil 
ounit= 22 

open (unit=ounit,file=outfil,status='new’,recl=idim2, 

+ iostat=ier,fomi='unforTnatted', access='sequential') 
if (ier.ne.O) then 

print *,'Error opening file ',outfil(25 : 35) 
stop 
endif 

c Open Landmask file 

mfile = drive//’:TOOLS:7/grdtyp//'MASK.DAT; 1 ’ 
print *, 'mask file = ', mfile 
munit = 18 

open (unit=munit,file=mfile, access='sequential', 

+ form=’unformatted',recl=idim2,status='old') 

c Get ice algorithm coefficients derived from tiepoints 

call usetie (c,ihem) 



c — start ice concentration algorithm 

do 1000 lp=l,idim3 

c Swap bytes in buffer file for Apple Macintosh 

do 800 inb=l,npar 
read(iunit) bbyte 

aternp(2)=btemp(l) 

atemp(l)=btemp(2) 

buff(inb)=abyte 

800 continue 

read(munit)(mask(inm),inm=l,idim2) 

c Average TBs into 25 km cells and determine ice concentrations 

call ice25 (buff, c,maxlen,idiml,idim2, ice, mask, icepar) 

c Write grid line of ice concentrations to file. 

write(ounit)(ice(iout),iout=l,idim2) 
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1000 continue 



c end ice algorithm 

close (unit=munit) 
close (unit=iunit) 
close (unit=ounit) 
close (unit=29) 

5000 continue 

6000 end 
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program locate 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c LOCATE - This program transforms I,J coordinates of an SSM/I grid cell 
c to latitude and longitude coordinates. This program provides 
c the inverse functions as well. LOCATE interfaces to the revised 
c forms of the subroutines, MAPXY and MAPLL. 
c 

c User-defined Paramters: 
c 

c gtype : Integer supplied by the user to describe one of the three 
c grid cell dimensions (12.5 km, 25.0 km, or 50.0 km), 

c 

c ihem : Integer supplied by the user to describe one of the two 
c polar regions (l=North , 2=South) 

c 

c itrans : Integer supplied by the user to describe the type of 
c transformation LOCATE will perform (l=I,J-to-Lat,Lon; 

c 2=Lat,Lon-to-I,J) 

c 

c i,j : Integers supplied by the user when itrans = 1. These 
c integers describe the position of a cell in an SSM/1 grid, 

c 

c lat.lon : Reals supplied by the user when itrans = 2. These 

c integers describe the latitude and longitude in an SSM/I 

c grid which LOCATE will transform to an I,J grid cell position, 

c Note: All latitudes and longitudes must be entered as 

c positive numbers! 

c 

c Internal: 
c 

c x,y 
c 
c 

c alat, 
c alon 
c 

c SGN 
c 
c 

c delta 
c 
c 
c 

c kk 
c 
c 

c SLAT : Standard latitude for the SSM/I grids is 70 degrees, 
c 

c numy : Number of lines in an SSM/I grid. This attribute varies 
c for each of the six grids. 



: Distance in kilometers from the origin of the grid 
(ie., pole). 



: Computed latitude and longitude returned from MAPXY. 

: Sign of the latitude (positive = north latitude, 
negative = south latitude) 

: Meridian offset for the SSM/I grids (0 degrees for 
the South Polar grids; 45 degrees for the North Polar 
grids. 

: Integer variable used for reorientation of the grid. The 
grid is 'flipped' in the Y direction for transformations. 
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c 

c cell : Size of the SSM/I grid ( 12.5 km, 25.0 km, 50.0 km) 
c 

c xydist : Distance from the origin of the grid in the cartesian plane, 
c The x-y coordinates for the edge of the lower left pixel 

c is (3850.0, 5350.0) for the northern grids and 

c (3950.0, 3950.0) for the southern grids, 

c 

c RE : Radius of the earth in kilometers, 
c 

c E : Eccentricity of the Hughes ellipsoid 
c 

c E2 : Eccentricity squared 

c 

c PI : Pi 

c Written by V.J.Troisi - January, 1990 

c Modified by W.A. Wright - July, 1991 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



real SLAT, E, RE, PI 

real alat,alon,x,y 

integer ihem 

real lat, Ion 

real SGN, delta 

integer numy(2,3) 

real cell(3), xydist(2,2) 

data numy / 896, 664, 448, 332, 224, 166 / 

data cell ' / 12.5 , 25.0, 50.0/ 

data xydist / 3850.0 , 5350.0 , 3950.0 , 3950 / 

SLAT = 70. 

RE = 6378.273 
E2 = .006693883 
PI = 3.141592654 
E = sqrt(E2) 
c 

c Query for the SSM/I grid cell size, 
c 

print *, ’Enter the grid cell dimension:’ 
print 1. 12.5 Km’ 
print *,’ 2. 25.0 Km’ 
print *,’3. 50.0 Km’ 
read *,gtype 
c 

c Query for polar region of interest, 
c 

print *, ’Enter the hemisphere of interest:’ 
print *,’ 1. North’ 
print *,’ 2. South' 
read *,ihem 
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c 

c Define the sign and meridian offset (delta) for the SSM/I grids, 
c 

if (ihem.eq.l) then 
SGN = 1.0 
delta = 45. 
else 

SGN =-1.0 
delta = 0.0 
endif 
c 

c Query for translation type, 
c 

print *, 'Enter one of the following transform functions:' 
print 1. Convert I,J to Latitude, Longitude' 
print 2. Convert Latitude, Longitude to I.J' 
read *, itrans 
c 

c Start translation 
c 

if (itrans.eq.l) then 
c 

c Obtain the I,J position of the grid cell to transform to Latitude 
c and Longitude 
c 

print *, 'Enter i,j:' 
read *,i,j 
c 

c Convert I,J pairs to x and y distances from origin. The grid will be 
c 'flipped' in the 'Y' direction, 
c 

x=((i- 1 )*cell(gtype))-(xydist( 1 ,ihem)-cell(gtype)/2.) 
kk=numy(ihem,gtype)-(j- 1 ) 

y=((kk-l)*cell(gtype))-(xydist(2,ihem)-cell(gtype)/2.) 

c 

c Transform x and y distances to latitude and longitude 
c 

call mapxy (x,y,alat,alon,SLAT,SGN,E,RE) 
c 

c Transform radians to degrees, 
c 

alon=alon* 1 80./PI 
alat=alat* 1 80./P1 
alon=alon-delta 
c 

c Convert longitude to positive degrees 
c 

if (alon.le.0.0) alon=alon+360. 
if (alon.ge. 360.0) alon=alon-360. 
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c Print the latitude and longitude for the center of the I,J cell, 
c 

print *,alat,alon 
pause 
else 
c 

c Obtain the latitude and longitude pair and transform to cell where 
c that pair is located, 
c 

print *, Enter latitude and longitude (positive values):' 
read *,lat,lon 
c 

c Transform degrees to radians 
c 

alat=abs(lat)*PI/l 80. 
alon=(lon+delta) *PI/ 1 80. 
c 

c Transform latitude and longitude to x and y distances from origin 
c 

call mapll (x,y,alat,alon,SLAT,SGN,E,RE) 
c 

c Convert x and y distances from origin to I,J pair (ii,jj) 
c 

ii=nint((x+xydist(l,ihem)-cell(gtype)/2.)/cell(gtype))+l 

jj=nint((y+xydist(2,ihem)-cell(gtype)/2.)/cell(gtype))+l 

c 

c Flip grid orientation in the 'Y' direction 
c 

kk=numy(ihem,gtype)-(jj- 1) 
c 

c Print the I,J location of the cell, 
c 

print *,ii,kk 
pause 

endif 

end 

SUBROUTINE MAPLL (X,Y,ALAT, ALONG, SLAT, SGN,E, RE) 

C$ 

C$ 

C$ DESCRIPTION: 

C$ 

C$ This subroutine converts from geodetic latitude and longitude to Polar 
C$ Stereographic (X,Y) coordinates for the polar regions. The equations 
C$ are from Snyder, J. P., 1982, Map Projections Used by the U.S. 

C$ Geological Survey, Geological Survey Bulletin 1532, U.S. Government 
C$ Printing Office. See JPL Technical Memorandum 3349-85- 101 for further 
C$ details. 

C$ 

C$ 
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I/O Description 



C$ ARGUMENTS: 

c$ 

C$ Variable Type 
C$ 

C$ ALAT REAL*4 I Geodetic Latitude (degrees, +90 to -90) 
C$ ALONG REAL*4 I Geodetic Longitude (degrees, 0 to 360) 
C$ X REALM O Polar Stereographic X Coordinate (km) 

C$ Y REALM O Polar Stereographic Y Coordinate (km) 

C$ 

C$ 

C$ 

C$ 

CS 
C$ 

C$ 



Written by C. S. Morris - April 29, 1985 
Revised by C. S. Morris - December 1 1, 1985 



Revised by V. J. Troisi - January 1990 
SGN - provides hemisphere dependency (+/- 1) 

REALM X,Y,ALAT,ALONG,E,E2,CDR,PI,SLAT,MC 

C S* ******************************************************************** 

c$ 

C$ DEFINITION OF CONSTANTS: 

C$ 

C$ Conversion constant from degrees to radians = 57.2957795 1 . 

CDR=57. 29577951 
E2=E*E 



C$ Pi=3. 141 592654. 

PI=3. 14 1592654 
C$ 

£<£********************************************************************* 

C Compute X and Y in grid coordinates. 

IF (ABS(ALAT).LT.Pl/2.) GOTO 250 
X=0.0 
Y=0.0 
GOTO 999 
250 CONTINUE 

T=TAN(PI/4.-ALAT/2.)/((l.-E*SIN(ALAT))/(l.+E*SIN(ALAT)))**(E/2.) 

IF (ABS(90.-SLAT).LT. 1 .E-5) THEN 
RHO=2.*RE*T/(( 1 .+E)**( 1 .+E)*( 1 .-E)**( 1 ,-E))**(E/2.) 

ELSE 

SL=SLAT*PI/1 80. 

TC=TAN(PI/4.-SL/2.)/(( l.-E*SIN(SL))/( 1 ,+E*SIN(SL)))**(E/2.) 
MC=COS(SL)/SQRT(1.0-E2*(SIN(SL)**2)) 

RHO=RE*MC*T/TC 
END IF 

Y=-RHO*SGN*COS(SGN*ALONG) 

X= RHO*SGN*SIN(SGN*ALONG) 

999 CONTINUE 
END 
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SUBROUTINE MAPXY (X,Y, ALAT, ALONG, SLA T,SGN,E, RE) 

C$ 

C$ 

C$ DESCRIPTION: 

CS 

C$ This subroutine converts from Polar Stereographic (X,Y) coordinates 
C$ to geodetic latitude and longitude for the polar regions. The equations 
C$ are from Snyder, J. P., 1982, Map Projections Used by the U.S. 

C$ Geological Survey, Geological Survey Bulletin 1532, U.S. Government 
C$ Printing Office. See JPL Technical Memorandum 3349-85- 101 for further 
C$ details. 

C$ 

c$ 

C$ ARGUMENTS: 

C$ 



Variable Type I/O Description 



X REAL*4 
Y REAL*4 
ALAT REAL*4 
ALONG REAL*4 



I Polar Stereographic X Coordinate (km) 

I Polar Stereographic Y Coordinate (km) 

O Geodetic Latitude (degrees, +90 to -90) 
O Geodetic Longitude (degrees, 0 to 360) 



C$ 

CS 
CS 
CS 
CS 
CS 
CS 
CS 
CS 
CS 

c$ 

CS 

c$ 

CS 

REAL*4X,Y, ALAT, ALONG, E,E2,CDR,PI 

(3S 

c$ 

CS DEFINITION OF CONSTANTS: 

c$ 

CS Conversion constant from degrees to radians = 57.29577951. 



Written by C. S. Morris - April 29, 1985 
Revised by C. S. Morris - December 1 1, 1985 

Revised by V. J. Troisi - January 1990 
SGN - provide hemisphere dependency (+/- 1) 



CDR=57. 2957795 1 
E2=E*E 

CS Pi=3. 141592654. 
PI=3. 141592654 



SL = SLAT*PI/180. 

200 RHO=SQRT(X**2+Y**2) 

IF (RHO.GT.0.1) GOTO 250 
ALAT=90.*SGN 
ALONG=0.0 
GOTO 999 

250 CM=COS(SL)/SQRT(1.0-E2*(SIN(SL)**2)) 
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T=TAN((PI/4.0)-(SL/(2. ())))/(( 1.0-E*SIN(SL))/ 
C(1.0+E*SIN(SL)))**(E/2.0) 

IF (ABS (SLAT-90.). LT.l.E-5) THEN 
T=RH0*SQRT((1.+E)**( 1 ,+E)*(l .-E)**(l .-E))/2./RE 
ELSE 

T=RHO*T/(RE*CM) 

END IF 

CHI=(PI/2.0)-2.0*ATAN(T) 

ALAT=CHl+((E2/2.0)+(5.0*E2**2.0/24.0)+(E2**3.0/12.0))*SIN(2*CHI)+ 

C((7.0*E2**2.0/48.0)+(29.0*E2**3/240.0))*SIN(4.0*CHI)+ 

C(7.0*E2**3.0/120.0)*SIN(6.0*CHI) 

AT AT=^nN*AT AT 

AL0NG=ATAN2(SGN*X,-SGN*Y) 

ALONG=SGN*ALONG 
999 CONTINUE 
END 



101 



program u, v calc 



* 

* This program calculates the u and v components of motion from a list of 

* user specified geodetic latitude and longitude postions. Results are output 

* as an ASCII data file in units of cm/s. Appropriate lines in code must be 

* modified to specify the filepath of input and output files. 

* 

* Variables: fname - input file name 

* infile - input file path 

* ofile - output file path 

* lat,lon - latitude and longitude expressed as real numbers 

* mlat - mean latitude between two positions 

* u, v - u and v components of motion 

* day - day of position (integer) 

* length - length in days of series 

* tmpd - time between postions (in days) 

* 



character* 10 fname 

character*70 ofile, infile 

real lat(60), lon(60), u(60), v(60),mlat 

integer day(60), length, tmpd 

c Construct input file path and open 

5 print *, 'Enter input file name' 
read (*,'(a)') fname 

print *, 'Enter number of days in series' 
read*, length 

print*, 'Enter separation time period (days)' 
read*, tmpd 

infile='Archs HardDisk:THESIS:Position Data:'//fname 
open (unit=15, file=infile,status='old',iostat=ier) 
if (ier.ne.O) then 

print*, infile, 'not accessible! <ret>' 
pause 
stop 
end if 

c Construct output file path and open 

ofile='Archs HardDisk:THES IS .Position Data:'//fname//'.uv' 
open (unit=18, file=ofile, status='new', form='formatted', 

+ iostat=ier) 
if (ier.ne.O)then 
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print*, 'Error opening ',ofile,' <ret>' 
pause 
stop 
endif 

c Read in data 

read (15,10) (day(i), lat(i), lon(i), i=l, length) 

c Compute u, v components (km/time period) 

do 100 i=l, length- 1 

v(i) = (lat(i+l) - lat(i)) * 111.6114 

mlat = (lat(i) + lat(i+l))/2. 

call gcircl(mlat,lon(i),mlat,lon(i+l),u(i)) 
if (lat(i+l) .It. lat(i)) u(i) = -u(i) 
print*, u(i), v(i) 

100 continue 
pause 

c Convert u, v components to cm/s 
do 200 i=l, length- 1 

u(i) = u(i) * ( 1 ./tmpd)* ( 1 000. ) * ( 1 00. ) * ( 1 ./24 .) * ( 1 ,/3 600. ) 
v(i) = v(i) * (l./tmpd)*(1000.)*(100.)*(l./24.)*(l./3600.) 

200 conunue 

c Write to file 

write (18,20) 'u', V 

write (18,20) 'em/s', 'em/s' 

write (18,30) (u(i), v(i), i=l, length-1) 

10 format (i6,3x,f5.2,2x,f6.2) 

20 format (lx, a,',',3x,a) 

30 format (lx, f6.2,',',3x,f6.2) 

end 
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subroutine gcircl(lat 1 ,lon 1 ,lat2,lon2,gcircl) 

c This program computes great circle distance. 

c Inputs are latitude/longitude pairs in decimal degrees. 

c Output is circumferential distance between the points in kilometers. 

implicit real *8 (a-z) 
real*8 dpr 

real*4 latl, lonl, lat2, lon2, gcircl 
integer again 

dpr = 57.29577951308232 



c Compute difference in longitudes 
t = (lon2 - lonl) / dpr 
c Compute r 

dsinr = dsin(t) * dcos(lat2/dpr) 
r = dasin(dsinr) 

c Compute k 

dsink = dsin(lat2 / dpr) / dcos(r) 
k = dasin(dsink) 

c Compute angle D 

dcosd = dcos(r) * dcos(k - (latl / dpr)) 
d = dacos(dcosd) * dpr 

c Compute distance based on 1 1 1.61 14 km/deg 
gcircl = d * 111.6114 
return 
end 
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Inputs 



c : Real array of coefficients. Refer to NODS User 
handbook for description of coefficients A-L. 

ihem : Integer denoting hemisphere (l=north, 2=south). Used 
in determining which set of coefficients to use. 



Internals 



subroutine usetie ( c , ihem ) 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c USETIE - This subroutine computes coefficients for the ice algorithm 
c using the NASA-approved tiepoints. 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



TB...W: Open water tiepoints. One for each of the three channels 
used in the ice concentration algorithm (19V, 19H, 37V) 

TB...F: Firstyear tiepoints. One for each of the three channels 
used in the ice concentration algorithm (19V, 19H, 37V) 

TB...M: Multiyear tiepoints. One for each of the three channels 
used in the ice concentration algorithm (19V, 19H, 37V) 



DW 

SW 



: Difference between component tiepoints for open water. 

: Sum of the component tiepoints for open water. 

DF : Difference between component tiepoints for first year ice. 
SF : Sum of the component tiepoints for first year ice. 

: Difference between component tiepoints for multiyear ice. 
: Sum of the component tiepoints for multiyear ice. 



DM 

SM 



Outputs 

c 



Array of coefficients which have been computed from 
tiepoints. 



real c(12) 
integer ihem 
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c Tiepoints for open water, first-year ice and multi-year ice 
c The first element of each tiepoint array represents the tiepoints 
c for the North Polar region; the second element is applied to the 
c South Polar region. 

c declare the two-element open water tiepoints 

real TB19VW(2), TB19HW(2), TB37VW(2) 

c declare the two-element first-year tiepoints 

real TB19VF(2), TB19HF(2), TB37VF(2) 

c declare the two-element multi-year tiepoints 

real TB19VM(2), TB19HM(2), TB37VM(2) 

real DW(2) , DF(2) , DM(2) 
real SW(2) , SF(2) , SM(2) 

c First component of each tiepoint pair is for Northern Hemisphere 
c Second component of each tiepoint pair is of Southern Hemisphere 
c These tiepoints approved for use by the NASA Sea Ice Algorithm 
c Working Group ( 1/90) 



data TB19VW / 175.3, 175.3/ 
data TB19HW / 97.7, 97.7/ 
data TB19VF/ 254.0, 251.2/ 
data TB19HF/ 236.0, 241.7/ 
data TB19VM / 223.2, 223.2/ 
data TB 19HM / 203.9, 203.9/ 
data TB37VW / 199.6, 199.6/ 
data TB37VF / 250.0, 248.3/ 
data TB37VM/ 186.3, 186.3/ 



DW(1) = TB19VW(ihem) - TB19HW(ihem) 
DF(1) = TB19VF(ihem) - TB19HF(ihem) 
DM(1) = TB19VM(ihem) - TB19HM(ihem) 
SW(1) = TB19VW(ihem) + TB19HW(ihem) 
SF(1) = TB19VF(ihem) + TB19HF(ihem) 
SM(1) = TB 19VM(ihem) + TB19HM(ihem) 
DW(2) = TB37VW(ihem) - TB19VW(ihem) 
DF(2) = TB37VF(ihem) - TB19VF(ihem) 
DM(2) = TB37VM(ihem) - TB 19VM(ihem) 
SW(2) = TB37VW(ihem) + TB19VW(ihem) 



106 



SF(2) = TB37VF(ihem) + TB19VF(ihem) 
SM(2) = TB37VM(ihem) + TB19VM(ihem) 



c Derive coefficients from tiepoints 

c(l) = DM(1)*DW(2) - DM(2)*DW(1) 

c(2) = DM(2)*SW(1) - DW(2)*SM(1) 

c(3) = DW(1)*SM(2) - DM(1)*SW(2) 

c(4) = SM(1)*SW(2) - SM(2)*SW(1) 

c(5) = DF(1)*(DM(2) - DW(2)) + DW(1)*(DF(2) - DM(2)) 

1+ DM(1)*(DW(2) - DF(2)) 

c(6) = DF(2)*(SM(1) - SW(1)) + DW(2)*(SF(1) - SM(1)) 
1+ DM(2)*(SW(1) - SF(1)) 

c(7) = DF(1)*(SW(2) - SM(2)) + DW(1)*(SM(2) - SF(2)) 
1+ DM(1)*(SF(2) - S\V(2)) 

c(8) = SF(2)*(SW(1) - SM(1)) + SW(2)*(SM(1) - SF(1)) 

1+ SM(2)*(SF(1)-SW(1)) 

c(9) = DF(2)*DW(1) - DF(1)*DW(2) 

c(10) = DW(2)*SF(1) - DF(2)*S\V(1) 

c(l 1) = SW(2)*DF(1) - DW(1)*SF(2) 

c(12) = SF(2)*S\V(1) - SF(1)*SW(2) 



return 

end 
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subroutine ice25 ( tb,c,maxlen,idiml,idim2,ice,mask,icepar ) 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c ICE25 - Compute ice concentrations from NSIDC CD-ROM data at 25 km 
c resolution, 

c 

c Software is modified from the Vince Troisi's original 50 km version by 
c W. A. Wright, Naval Postgraduate School, 

c 

c Algorithm 

c 

c The ice concentration algorithm is an adaptation of the 

c algorithm designed by the NODS Software Group (refer to the 

c NODS User Handbook, Appendix CC, pp 22-26). The algorithm 
c uses the latest SSM/I tiepoints which were determined by Don 
c Cavalieri and Koni Steffen, 

c 

c A weather filter (GR > .05) is applied to correct for false 
c concentrations computed in regions of open water, 
c 

c Computed concentrations greater than 120% and less than -20% 
c are considered suspect and the corresponding cell is flagged 
c as missing (-99). 

c 

c Inputs 
c 

c tb : buffer of TB cell values. This buffer includes all 
c SSM/I TB parameters excluding the 85GHz channels, 

c 

c idiml : The number of parameters in a SSM/I cell (ie., 5) 
c 

c idim2 : The number of cells in a SSM/I Tb grid line, 
c 

c maxlen: The maximum length of tb 
c 

c mask : A line from the landmask for 25 km grid. The mask 
c line corresponds to the ice grid line to be 

c computed. (l=land; 0=water) 

c 

c icepar: The ice concentration type to compute (l=total, 
c 2=multiyear) 

c 
c 

c Internals 
c 

c c : Array of coefficients used to compute ice concentrations, 
c Refer to NODS Handbook Volue II. Appendix CC pp 23-24. 

c 

c ioff : Offset pointers used to point to Tb values used in 
c this ice algorithm (19GHz V, 19GHz H, 37GHz V). 

c 

c ipos : position in tb array which is to be summed 
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c 

c sum : Real register used for summing four adjacent Tb values 
c 

c count : The number of Tb values not flagged as missing. A Tb value 
c of zero (0) is a missing Tb value, 

c 

c fyr : Integer register which stores the calculated first-year 
c ice concentration 

c 



myr : Integer register which stores the calculated multi-year 
ice concentration 



: Integer register which stores the calculated total ice 
concentration 



fyrb : Two - byte buffer equivalenced to fyr register. Only the 
first byte is the significant part of the register. 

myrb : Two - byte buffer equivalenced to myr register. Only the 
first byte is the significant part of the register. 



c 
c 
c 

c tot 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c Outputs 
c 

c ice 
c 
c 
c 
c 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



GR 



PR 



: Gradient ratio 
: Polarization ratio 



: Line of 25 km ice concentration grid. Values are 
represented as a single byte percentage ranging 
from 0 - 100%. Cells over land flagged with -88; 
cells for which concentrations were suspect or 
not computed are flagged with -99. 



integer idiml,idim2 
integer maxlen 
integer*2 tb(maxlen*idiml) 



c Offset from start of cell to three parameters used in ice 
c calculations for two lines of data 



integer*2 ioff(3) 

integer*2 fyr,myr,tot 

integer* 1 fyrb(2),myrb(2) 

integer* 1 ice(idim2),mask(idim2),fyrice,myrice 

logical GOOD 

c Three TB parameters for ice concentrations 
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real BT19V,BT19H,BT37V 

real c(12) 

real GR,PR,PRGR 

real ANF,ANM,DD 

equivalence (fyr, fyrb(l)) 
equivalence (myr, myrb(l)) 

c Offset pointers for Brightness temperature cells used in algorithm 
data ioff /0, 1 ,3/ 
c Begin processing buffer 
ii = 0 

do 1000 i= l,idiml*idim2,idiml 
c Increment pointer for next Ice cell. 
ii=ii+l 

c Calculate Ice Concentradon 
GOOD = .true. 

c Set ice pixel to land(-88) if over land 

if (mask(ii).eq.l) then 
ice(ii)=-88 
else 

c Determine Ice Concentration of the ice pixel 

do 80 1=1,3 

ipos=i+ioff(l) 

tbval = (float(tb(ipos)))/10. 

if(l.eq.l) BT19V = tbval 
if(l.eq.2) BT19H = tbval 
if(l.eq.3) BT37V = tbval 
80 continue 
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c Set ice cell to missing if any TB parameters are missing 

if (BT19V.EQ.0.0 .or. BT19H.EQ.0.0 .or. 

* BT37V.EQ.0.0 ) then 
fyrice = -99 
myrice = -99 
ice(ii) = -99 

else 

c Compute open water weather filter 

GR = (BT37V - BT19V)/(BT37V + BT19V) 

c 

c Set ice cell to 0 % concentration if w'eather filter > .05 
c Filter value of 0.05 may be altered at a future date! 

c 

if (GR ,gt. .05) then 
fyrice = 0 
myrice = 0 
ice(ii) = 0 
else 

c Apply noniterative ice algorithm 

c Compute polarization ration PR 

PR = (BT19V - BT19H)/(BT19V + BT19H) 

PRGR = PR*GR 

ANF = c(l) + c(2)*PR + c(3)*GR + c(4)*PRGR 
ANM = c(9) + c(10)*PR + c( 1 1)*GR + c(12)*PRGR 
DD = c(5) + c(6)*PR + c(7)*GR + c(8)*PRGR 

c Compute First-Year, Multi-Year, and Total Ice Concentrations 

fyr = ifix(100.0*(ANF/DD)) 
myr = ifix(100.0*(ANM/DD)) 
tot = fyr + myr 

c 

c Filter out values outside the range of -20 to 120 
c This filter may be altered at a future date 

c 

if (fyr .It. -20 .or. 

* myr .It. -20 .or. 

* tot .It. -20 .or. 

* fyr .gt. 120 .or. 

* myr .gt. 120 .or. 

* tot .gt. 120) then 
fyrice = -99 
myrice = -99 
ice(ii) = -99 
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GOOD = .false, 
endif 

c Place concentration values into byte data representations 

if (GOOD) then 
fyrice = fyrb(2) 
myrice = myrb(2) 

if (icepar.eq.l) ice(ii) = fyrice + myrice 
if (icepar.eq.2) ice(ii) = myrice 

c 

c Set concentration to 0 if computed concentration between -20 and 0 
c Set concentration to 100 if computed concentration between 100 and 120 

c 

if (ice(ii) .It. 0) ice(ii) = 0 
if (ice(ii) .gt. 100) ice(ii) = 100 
endif 

GOOD =.true. 
endif 
endif 
endif 

1000 continue 

c Return computed ice concentration array, ice 

return 

end 
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